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A  detailed  derivation  of  improved  ocean  tidal  equations  in 
continuous  (COTEs)  and  discrete  (DOTEs)  forms  is  pre¬ 
sented.  These  equations  feature  the  Boussinesq  linear  eddy 
dissipation  tew  with  a  novel  eddy  viscosity  that  depends  on  die 
lateral  mesh  area,  i jt.,  on  mesh  size  and  ocean  depth.  Analo¬ 
gously,  the  linear  law  of  bottom  friction  Is  used  with  a  new 
bottom  friction  coefficient  depending  on  the  bottom  mesh  area. 
The  primary  astronomical  tide-generating  potential  is  modified 
by  secondary  effects  due  to  the  oceanic  and  terrestrial  tides.  The 
fully  linearized  equations  are  defined  in  a  single-layer  ocean 
basin  of  realistic  bathymetry  varying  from  SO  m  to  7,000  m. 
The  DOTEs  are  set  up  on  a  1°  by  1*  spherically  graded  grid 
system,  using  central  finite  differences  in  connection  with  Rich¬ 
ardson’s  staggered  computation  scheme.  Mixed  single-step  finite 
differences  in  time  are  introduced,  which  enhance  decay,  disper¬ 
sion,  and  stability  properties  of  the  DOTEs  and  facilitate — in 
Part  II  of  this  paper — a  unique  hydrodynamical  interpolation 
of  empirical  tide  data.  The  purely  hydrodynamical  modeling  is 
completed  by  imposing  boundary  conditions  consisting  of  no¬ 
flow  across  and  free-slip  along  the  mathematical  ocean  shore¬ 
lines.  Shortcomings  of  the  constructed  preliminary  M2  ocean 
tide  charts  are  briefly  discussed.  Needed  improvements  of  the 
model  are  left  to  Part  II. 
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Up  to  recent  years  practical  interest  in  ocean  tides  was  essentially 
confined  to  coastal  waters.  With  the  advancement  of  science  and 
technology,  the  need  for  extremely  accurate  tide  predictions  in  all 
the  world  oceans  has  become  an  urgent  problem  in  applications  to 
geophysics,  geodesy,  oceanography,  meteorology,  astronomy,  and 
space  technology.  In  the  past  two  decades  considerable  progress  in 
the  qualitative  mapping  of  global  ocean  tides  has  been  made  by 
hydrodynamical  and  empirical  methods.  Yet,  the  constructed  tidal 
charts  vary  considerably  over  large  ocean  areas  from  investigator 
to  investigator,  and  tide  predictions  fall  considerably  short  of  the 
desired  accuracy.  A  review  of  the  numerical  work  was  recently 
published  by  Hendershott  (1977).  A  comprehensive  review  of  iiic 
major  highlights  in  ocean  tidal  research  from  the  beginning  to  the 
latest  results  is  in  preparation  by  the  author  (Schwiderski,  1980a). 

In  the  present  paper  the  well-known  hydrodynamical-numeri- 
cal  method  developed  by  Hansen  ( 1966)  and  analyzed  and  tested 
by  Zahel  (1970,  1973,  1975)  and  Estes  (1975,  1977)  will  be 
improved  with  respect  to  eddy  dissipation,  bottom  friction,  grid 
system,  and  finite  differencing  technique.  The  newly  derived  dis¬ 
crete  ocean  tidal  equations  (DOTEs)  will  be  analyzed  and  applied 
to  construct  a  preliminary  Ms-tide  model,  using  the  ocean  bathyme¬ 
try  data  collected  by  Smith  et  al.  (1966).  The  results  of  this  purely 
hydrodynamical  technique  will  be  critically  evaluated  and  needed 
improvements  indicated.  In  a  subsequent  paper  (Part  II,  Schwider¬ 
ski,  1979b)  tee  shortcomings  of  the  present  theoretical  model  will 
be  eliminated  by  using  a  hydrodynamically  defined  bathymetry  and 
a  novel  hydrodynamical  interpolation  of  empirical  tidal  constants 
collected  around  tee  world. 


Derivation  of  Continuous  Ocean  Tide  Equations 
The  Navier-Stokes  Equations  of  Averaged  Turbulent  Flow 
Because  of  tee  enormous  dimensions  of  tee  world  oceans,  their 
hydrodynamical  tidal  motion,  generated  by  the  attraction  of  the 
moon  and/or  sun,  must  be  considered  entirely  turbulent.  Accord¬ 
ingly,  any  comprehensive  modeling  of  oceanic  tidal  currents  should 
begin  with  tee  complete  Navier-Stokes  equations  of  averaged  tur- 
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buleat  motions  of  a  viscous,  incompressible  fluid  including  all 
unknown  Reynolds  stresses  (see,  e.g.,  Schlichting,  1968),  which 
must  provide  the  major  friction  stress  to  keep  the  intrinsic  super¬ 
critical  instability  of  the  flow  under  control.  In  rotating  spherical 
polar  coordinates  (A,  4,  r),  these  equations  may  be  written  in  the 
form  (see,  e.g.,  Whitaker,  1968) : 

«,  +  —2—rUK  +  Vtt  +  l£V_i£Vtan^+i£HL 

r-cos  d»  r  r  r  r 
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III  these  equations,  all  subscripts  denote  indicated  partial  deriva¬ 
tives,  while  all  superscripts  denote  indicated  components  of  the 
turbulent  dissipation  vector  /  or  the  Reynolds  stress  tensor  r.  Fur¬ 
thermore,  the  following  notations  are  used: 

/  =  universal  time  in  sec 

*  =  east  longitude 

0  =  |  -  9  =  north  latitude  (0  =  colatitude) 

r  —  polar  radius  in  m 

(u,  v,  w)=  Coot,  lMJilii,  and  radial  velocities,  respectively,  in 
m/sec  (averaged) 
p  —  pressure  t averaged) 

q  =  total  tide-generating  potential 

/  =  dissipation  vector 

t  =  Reynolds  stress  tensor  (to  be  specified) 

o  =  0.72722  x  10"4  sec"1  =  earth  angular  velocity 
p  ~  103  kg/m3  =  density  of  sea  water 

G  =9.81  m/sec*  =  gravity  acceleration 

The  equations  of  turbulent  “mean"  motion  are  obtained  by  a 
formal  time-averaging  procedure  applied  to  die  Navier-Stokes 
equations  of  viscous  laminar  flow  that  leads  to  the  concepts  of 
“averaged"  velocities  and  pressure,  as  well  as  to  the  unknown 
Reynolds  stress  tensor,  r,  containing  the  filtered  out,  fluctuating 
velocity  residuals  in  quadratic  form  (see,  e.g.,  Schlichting,  1968). 
A  discussion  of  the  physical  meaning  of  these  and  following  turbu¬ 
lence  motions  may  be  postponed  to  the  section,  “Discrete  Versus 
Continuous  Ocean-Tide  Equations."  The  momentum  equations 
(Equations  1)  maintain  the  momentum  balance  between  tfie  famil¬ 
iar  kinetic  forces  on  the  left  side  and  the  forces  of  potential  (tide, 
gravity,  and  pressure),  Coriolis,  centrifugal  acceleration,  and  dis¬ 
sipation  on  the  right  side.  The  continuity  equation  (Equation  2) 
expresses  die  condition  of  conservation  of  incompressible  mass. 
The  Reynolds  stress  tensor,  r,  and  the  total  tide-generating  poten¬ 
tial,  q,  will  be  specified  in  the  next  two  sections.  The  global  oceans 
at  hydrostatic  rest  may  be  described  by  the  following  surface  and 
bottom  boundary  conditions: 
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( 1 )  r  =  R  =  0.637  x  10T  m  =  spherical  (geoidal  rest)  sea  sur¬ 

face  (all  geoidal  undulations  are 
neglected  without  any  significant 
loss  of  accuracy), 

(2 )  r  =  R  —  H(a,  4)  =  sea-bottom  relief,  where 

(3)  H  =  H(x,  4>)  =  realistic  ocean  depth  in  m  H(A,  4)  =  0 

for  land  points,  see  section  “The  1 c  by  1° 
Graded  Grid  System  and  Bathymetry,” 
and 

(4)  p  =  P  -  Gpz  =  hydrostatic  sea-pressure  distribution  with 

constant  (arbitrary)  sea-surface  pressure 
P,  where 

(5) z  =  r-  R=  new  depth  variable,  so  that  z  =  0  denotes 

r^R  (see  Figure  1). 

Due  to  the  time-dependent,  tide-generating  potential,  q,  acting 
on  the  ocean  and  solid  earth,  die  hydrostatic  conditions  (l)-(S) 
are  altered  to  the  following  hydrodynamical  boundary  values  (see 
Figure  1): 

2  (r) 

i 

2*»  £*  •  TIDAL  SURFACE 
Z  •  0  ■  GEOIDAL  SURFACE 

}  H  •  DEPTH 

Zb'(b-H 'BOTTOM  TIDAL  RELIEF 
Z'-H'  BOTTOM  REST  RELIEF 

1.  Baith-odean  tidal  tateructfon:  { «  ocean  tide,  f*  =  earth  tide,  {•  m 
•urface  tide,  I*  *  bottom  tide,  ss  f*  —  f*  =  earth  dip-response  to  ocean 
tide,  and  H  ss  ocean  depth. 
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(a)  z*  =  C  =  C  (K  0,  t)  —  total  sea-surface  tidal  elevation 
over  geoidz  =  0; 

(b)  zb  =  (X,  p,  t)  -  H(A.  0)  =  sea-bottom  tidal  relief,  where 

(c)  (A,  0,  t)  =  total  bottom  tidal  elevation  over  rest 
relief  z  =  -H(A,0); 

(d)  {  =  (  (K  0,  t)  =  C  ~  =  ocean  tidal  elevation  (measured 

by  bottom  tidal  pressure  gauges)  to  be  modeled; 

(e)  C  =  C  (X,  0,  t)  =  earth  tidal  elevation  to  be  specified 
(Equations  9) ; 

(f)  £*°  =  (A,  0,  t)=  £*  -  C*  =  earth  dip-response  to  oceanic 

tidal  load,  ( (Equation  10); 

(g)  p*  =  P  =  constant  surface  pressure  (no  atmospheric  pres¬ 
sure  considered) ; 

(h)  w*  =  £  -  surface  radial  velocity; 

(i)  t*  =  o  =  no  surface  stress  (frce-slip,  no  wind  force  con¬ 
sidered); 

(j)  (ub,  v",  w")  •  V(H  +  z)~Cl  =  no  flow  across  ocean  bot¬ 
tom  (V  =  gradient  vector) ; 

(k)  r*  =  *B(ub,  v*)  =  bottom  stress  vector  specified  by  invok¬ 
ing  the  linear  taw  of  bottom  friction  with  the  coefficient 

B  —  bpcos  0 = FL*  a*cos  0  (4a) 

depending  on  the  cell  bottom  area  L’pcos  4,  where 

(l)  L  =  chosen  equatorial  mesh  size  (section  The  1°  x  1° 
Graded  Grid  System  and  Bathymetry”), 

(m)  m  =  mesh  “grading”  parameter  (section  “The  1°  x  1° 
Graded  Grid  System  and  Bathymetry”),  and 

(n)  b  —  EL2  =  uniform  bottom  friction  parameter,  which 

must  be  determined  by  trial-and -error  computations  for 
best  results.  The  final  M2  tide  (Schwiderski,  1979b)  was 
computed  with  b  =  0.01  m/sec.  (4b) 

(o)  Lateral  boundary  data  will  be  specified  in  die  section, 
“Lateral-Boundary,  Initial,  and  Final  Data.” 


In  contrast  to  Zahel  (1970,  1973,  1975)  and  Estes  (1975, 
1977),  who  used  the  nonlinear  law  of  bottom  friction,  the  linear 
law  is  hare  preferred  as  it  is  more  consistent  with  die  overall  Stokes 
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skrw-motioo  assumptions  (see,  e.g.,  Schlichting,  1968)  character¬ 
izing  the  present  tide  model  (section,  “The  Continuous  Ocean-Tide 
Equations  (COTEs).”  Especially  the  computations  by  Estes  (1975, 
Figure  6)  indicate  no  need  for  any  nonlinear  friction  terms,  at  least 
not  in  deep  ocean  areas.  Pekeris  and  Accad  (1969)  used  the  same 
linear  law  of  bottom  friction  and  called  attention  to  earlier  work 
by  Grace  (1931)  in  the  Gulf  of  Suez.  Grace  applied  both  laws  and 
experienced  a  slight  preference  for  the  linear  law. 

In  deviation  from  Pekeris  and  Accad,  the  present  bottom  fric¬ 
tion  coefficient  is  assumed  independent  of  ocean  depth.  Since  those 
authors  placed  their  artificial  boundaries  of  the  world  oceans  at 
the  1,000-m  depth  level,  it  was  plausible  to  restrict  any  significant 
bottom  friction  to  those  boundaries.  This  was  accomplished  by 
simply  assuming  a  bottom  friction  coefficient  B  that  depended  in¬ 
versely  on  the  ocean  depth.  In  the  present  tide  model,  which  uses 
(section,  “The  l°xl°  Graded  Grid  System  and  Bathymetry”),  a 
bathymetry  that  varies  from  50  m  to  7,000  m  (or  10  m  to  7,000  m 
in  Schwiderski,  1979b),  this  assumption  was  found  unjustified. 
Bottom  friction  coefficients  of  the  form  B~H*  have  been  tested  for 
e  =  -1  (Pekeris  and  Accad  assumption),  e  —  -V4,  and  e  -  o.  The 
hater  case  appeared  to  yield  best  results. 

The  assumed  dependence  of  the  bottom  friction  coefficient  B  on 
the  bottom  mesh-cell  area  (Equation  4a)  is  in  agreement  with  the 
novel  law  of  eddy  viscosity,  which  is  introduced  below  in  die  next 
section  and  which  is  further  discussed  in  the  section,  “Discrete 
Versus  Continuous  Ocean-Tide  Equations.”  In  fact,  as  can  be  seen 
from  Equations  28,  the  bottom  friction  coefficient  is  essentially  de¬ 
termined  by  the  vertical  eddy  viscosity.  It  may  be  noted  that  tire 
numerical  values  of  the  bottom  friction  coefficient  (Equations  4a 
and  b)  used  in  the  present  model  fall  considerably  below  those  used 
by  Pekeris  and  Accad  in  their  corresponding  domains,  where  fric¬ 
tion  is  essentially  restricted  to  the  boundary  line. 

Brynniifii  5  frames  mrf  Flrfrfy  rHas^w  ninn 
The  formally  derived  unknown  Reynolds  stress  tensor,  r  (preced¬ 
ing  section),  may  be  specified  by  the  original  Boossinesq  (1877) 
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substitution  (see  section,  “The  l°x  1°  Graded  Grid  System  and 
Bathymetry”);  that  is,  by  the  laminar,  viscous  stress  tensor  (see, 
e.g.,  Schlichting,  1968;  Whitaker,  1968): 


'“  =  2"4  7*“*+r]* 

(5a) 

w[K+t]  • 

(5b) 

r"  =  2pAwr , 

(5c) 

rcos*  ( u  \  1  vi 

p  L  r  Vcos  '  r  cos  4  XJ  ’ 

(5d) 

r»  =  r»  =  fA[r  (f)r+fcos,n]. 

(5e) 

* 

II 

I 

II 

X 

n  i*- 

(5f) 

In  this  substitution,  the  ordinary  kinetic  molecular  viscosity  is  re¬ 
placed  by  the  so-called  “eddy  viscosity”  A  (momentum  austausch, 
exchange,  or  mixing  coefficient),  which  remains  to  be  modeled  to 
represent  the  true  turbulent  flow  characteristics  at  hand  as  closely 
as  possible. 

In  order  to  achieve  a  greater  modeling  flexibility,  it  is  customary 
to  divide  flte  eddy  viscosity  into  a  ‘Vertical”  eddy  viscosity  associ¬ 
ated  with  vertical  shear  and  into  a  “horizontal”  (lateral)  eddy 
viscosity  associated  with  horizontal  shear.  Since  east  and  north 
dimensions  are  usually  nearly  equal,  no  need  for  two  lateral  eddy 
viscosities  has  ever  been  encountered.  In  single-layer  ocean  mod¬ 
els  like  tbe  present  one  (section,  “The  Continuous  Ocean-Tide 
Equations  (COTEs),”  the  separate  treatment  of  the  vertical  eddy 
viscosity  needs  no  explicit  specification,  because  (Equations  28) 
the  viscosity  becomes  an  important  part  of  the  bottom-friction  co¬ 
efficient,  B,  defined  by  Equations  4a  and  4b. 

Based  on  the  author’s  extensive  computer  experiments  and  on 
the  physical  arguments  elaborated  on  in  the  sections,  “Stability 
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Analysis”  and  “Discrete  Versus  Continuous  Ocean-Tide  Equa¬ 
tions,”  the  (horizontal)  eddy  viscosity.  A,  may  now  be  specified  by 

A  —  yLf/(A,  4>)  ( 1  +  mcos0).  (6a) 

In  this  definition,  L  and  m  denote  the  mesh  size  and  grading  param¬ 
eters  introduced  in  the  preceding  section  (l,m).  Accordingly,  the 
new  eddy  viscosity  depends  on  the  mean  lateral  cross-section  area 
H  x  L(1  +  m  cos  <t>)/2  of  the  flow  cell  of  depth  H  and  average 
north-east  mesh  size  L(1  +  m  cos  The  remaining  reduced 

eddy-dissipation  coefficient,  a  (in  sec),  must  be  subjected  to  trial- 
and-error  computations  in  order  to  achieve  best  results  uniformly 
over  all  oceans.  (Different  a-values  for  the  Pacific,  North  and  South 
Atlantic,  and  Indian  oceans  were  also  tested  without  significant 
effects. ) 

At  this  point,  it  may  be  mentioned  that,  following  Zahel  (1970, 
1973),  the  author  used  in  an  exploratory  tide  model  a  constant 
eddy  viscosity  uniformly  over  all  oceans.  The  results  were  rather 
disappointing,  especially  for  shallow  ocean-shelf  areas,  where  in¬ 
tolerably  low  tidal  amplitudes  were  computed.  Subsequent  com¬ 
puter  experiments  with  eddy  viscosities  of  the  form  A~H*  with 
e  —  o,  1/2,  1,  and  3/2  indicated  best  results  for  e  =  1,  as  in  Equa¬ 
tion  (6a).  Some  dependence  of  the  eddy  viscosity  on  the  mesh  size 
considered  was  earlier  noticed  by  Cox  (1970),  Friedrich  (1970), 
Holland  and  Hirschman  (1972),  and  Zahel  (1975).  Indeed,  eddy 
viscosities  ranging  from  A  =  10  to  10,1cm2/scc  have  been  reported 
as  required  by  many  researchers  investigating  different  problems. 
Obviously,  those  huge  variations  of  the  eddy  viscosity  can  be  ex¬ 
plained  by  Equation  (6a)  by  the  strong  variations  of  the  mesh  size 
L  and  the  depth  H  considered  by  different  analysts.  In  the  present 
case  with  H  varying  from  10  m  to  7,000  m  (section,  “The  l°xr 
Graded  Grid  System  and  Bathymetry”)  the  eddy  viscosity  was 
found  to  vary  between  (section,  “Stability  Analysis”) 

1.3*  10* ~-<  A  <  1.3  •  10“—, 


(6b) 
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which  fits  the  customary  range  very  well.  It  is  perhaps  fortuitous, 
yet  interesting,  to  note  that  for  a  cubic  cell  of  about  1-in.  mesh  size 
the  eddy  viscosity,  A  (Equations  6a,  b),  reduces  to  the  value  of  the 
molecular  viscosity  of  water. 

The  Total  Tide-Generating  Potential 

The  total  tide-producing  potential,  q  (section,  “The  Navier-Stokes 
Equations  of  Averaged  Turbulent  Flow”),  may  be  expressed  in  the 
form 


q=-G(v  +  n),  (7) 

where  Gv  is  the  “primary”  astronomical  potential  directly  propor¬ 
tional  to  Newton’s  equilibrium  tide  7.  The  remaining  “secondary” 
potential  Gv'  can  be  expanded  into  its  three  major  parts 

Gv'  =  G(v°+  v'-v**),  (8) 

which  reveal  their  corresponding  origins  (Figure  1,  section,  “The 
Navier-Stokes  Equations  of  Averaged  Turbulent  Flow”  (d),  (e), 
(f) ),  respectively: 

Gv0  =  gravity  potential  perturbation  due  to  the  ocean  tide  (, 
Gi j*  =  gravity  potential  perturbation  due  to  the  earth,  tide  ?, 
Gv*°  —  gravity  potential  perturbation  due  to  the  earth  dip-response 
C*  caused  by  die  ocean  load  ((-tide). 

Attention  to  the  significance  of  the  earth-ocean  tidal  interac¬ 
tions  manifested  in  the  five  quantities  {*,  v0,  F®,  v”,  and  v°  was 
ceiled  for  by  Farrell  (1972a,  b;  1973).  In  ocean-tide  models,  it  is 
sufficient  to  use  the  following  approximate  relations: 


and 


(*  **  0.61?,  v*  “*  0.30i», 

(9a) 

r-<~0.31n. 

(9b) 

tr  -  +  v*  ~  o.ioc 

(10) 

Ocean  Tides,  Part  I:  Global  Ocean  Tidal  Equations 


171 


The  first  relations  (Equations  9)  used  by  Hendershott  (1972, 
1975),  Zahel  (1975),  Estes  (1977),  and  others,  can  be  justified 
by  second-order  Love  number  approximations.  The  second  relation 
(Equation  10)  has  been  suggested  by  Pekeris  (presentation' given 
in  1977  at  the  U.S.  Office  of  Naval  Research),  who  also  recom¬ 
mended  the  factor  0.10  after  evaluating  the  Green’s  function  repre¬ 
sentation  of  the  three  oceanic  tidal  load  effects  (C°,  r;M,  *j°)  derived 
by  Farrell.  Apparently,  the  suggestion  by  Pekeris  (Equation  10)  is 
physically  just  as  plausible  as  the  accepted  approximation  (Equa¬ 
tions  9).  In  Equations  9a  and  b,  it  is  assumed  that  the  solid-earth 
tide,  r,  and  its  subsequent  gravity  perturbation,  Gye,  are  essentially 
instantaneous  responses  to  the  tide-generating  potential,  q.  Analo¬ 
gously,  in  Equation  10,  it  is  assumed  that  the  solid-earth  dip,  F°, 
and  the  gravity  perturbations,  and  v°,  are  almost  instantaneous 
responses  to  the  ocean’s  tidal  load,  t.  It  may  be  mentioned  that  the 
author  conducted  extensive  computer  experiments  using  die  factors 
0.00,  0.08,  and  0.12  instead  of  0.10  in  Equation  10.  The  last  two 
factors  produced  no  noteworthy  alterations  of  the  results.  The  first 
factor  (0.00)  obviously  deletes  all  oceanic  tidal-load  effects  to 
which  critically  large  significance  has  been  attached  by  Hender¬ 
shott  (1972,  1975).  The  author’s  computations  supported  the 
marginal  effects  of  oceanic  tidal  loading  found  by  Estes  (1977). 

Following  Thomson  (1868),  Darwin  (1883),  Doodson  (1921), 
and  Cartwright  and  Taylor  (1971 ),  the  primary  astronomical  tide¬ 
generating  potential,  G>?,  or,  equivalently,  the  equilibrium  tide,  v, 
may  be  expanded  into  a  series  (see,  e.g.,  Dietrich,  1963)  of  “har¬ 
monic  components"  (constituents),  n  with  a  nonharmonic  fre¬ 
quency  spectrum 

v  -  ^-rr-Q  <*>,  t).  (11) 

In  the  final  analysis  of  ocean  tides,  the  following  three  major  spe 
cies  will  have  to  be  considered : 

(a)  Semidiurnal  Equilibrium  Tides 

v  =  2:72  =  Kcos2  <f>  cos(<rt +  2a  +  x)  (12) 
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(b)  Diurnal  Equilibrium  Tides 

v  =  l:^i  =  Asm  2<f>  cos(vt  +  A  +  x)  (13) 

(c)  Long-Period  Equilibrium  Tides 

v  =  0:i?o  =  K(3  cos2  <t>  —  2)  cos(of  +  x)  (14) 

In  Equations  12-14,  the  constants  (K,  a,  x)  denote  K  =  amplitude 
of  equilibrium  tide  (in  m),  a  =  frequency  of  equilibrium  tide 
(in  sec-1),  x  =  astronomical  argument  of  equilibrium  tide  (in  rad). 
In  Table  1,  these  constants  are  listed  for  the  major  tidal  modes 


Table  1 

Constants  of  ma|or  tidal  modes. 


Tidal  Mode 

K«(m)  <rb(10-‘/sec)  x'Tdeg) 

Semidiurnal  Species 

M,  =  Principal  Lunar 

0.242  334 

1.405  19  2h0“  -  2s0 

S3  —  Principal  Solar 

0.112  841 

1.454  44  0 

N3  =  Elliptical  Lunar 

0.046  398 

1.378  80  2h0-3s0  +  p0 

K3  =  Declination  Luni-Solar  0.030  704 

1.458  42  2h0 

Diurnal  Species 

K,  =  Declination  Luni-Solar  0.141  565 

0.729  21  h0  +  90 

O,  =  Principal  Lunar 

0.100  514 

0.675  98  h0  —  2s0  —  90 

P,  =  Principal  Solar 

0.046  843 

0.725  23  -h0 -90 

Q,  =  Elliptical  Lunar 

0.019  256 

0.649  59  h0  -  2s0  +  Po  -  90 

Long-Period  Species 

Mf  =  Fortnightly  Lunar 

0.041  742 

0.053  234  2s0 

Mm  =  Monthly  Lunar 

0.022  026 

0.026  392  s0  —  Po 

Ssa  =  Semiannual  Solar 

0.019  446 

0.003  982  2h0 

a  K  =  amplitude  of  the  partial  tide. 

&  a  =  frequency  of  the  partial  tide, 
c  x  —  astronomical  argument  of  the  partial  tide. 

4  (ho.  Sg,  Pq)  =  mean  longitudes  of  sun,  moon,  and  lunar  perigee  at  Greenwich  midnight; 
h0  =  279.006  68+  36  000.  768  930  *85  T  +  3.08  •  10-«T* 

So  =  270.434  358  +  481  287.883  141  37T  -  0.001  1331s  +  1.9  •  10-«T*. 

Po  =  334.329  853  +  4  089.034  032  957  5T  -  0.010  325T*  -  1.2  •  10-*n, 

where 

T  =  (27  392.500  528  +  1.000  000  035  6DJ/38  525. 

D  =  d  +  365  (y  -  1975)  +  Int  ((y  -  1975)/4J. 
d  sb  day  number  of  year  (d  =  1  for  January  1), 
y  >  1975  =  year  number, 
and 

Int  (x]  s  integral  pert  of  x. 
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with  amplitudes  larger  than  4%  of  the  leading  semidiurnal  moon 
(Ms)  tide.  The  daily  astronomical  argument,  x»  can  be  neglected 
in  the  following  construction  of  the  oceanic  tidal  modes 

t  -  £(a,  4)  cos(o/  +  x  ~  S(x,4))  (15a) 

corresponding  to  the  considered  mode  of  the  equilibrium  tide,  v  = 
Vv  According  to  Equation  15a,  only  the  “harmonic  constants” 

£  =  £(A,  4>)  =  tidal  amplitude  (in  m) 
and  (15b) 

S  =  3(a,  <t>)  =  Greenwich  phase  (in  rad) 

need  to  be  found,  provided  the  tide  model  is  linear  or  almost  linear. 

The  Continuous  Ocean-Tide  Equations  ( COTEs) 

In  addition  to  the  simplifying  assumptions  made  in  the  three  pre¬ 
ceding  sections  concerning  bottom  friction  ((k)-(n)),  eddy  dissi¬ 
pation,  and  total  tide-generating  potential,  the  following  simplifica¬ 
tions  may  be  invoked: 

(a)  Hydrostatic  pressure  assumption  in  Equations  1 ;  i.e., 

( 1 )  Neglect  all  quadratic  inertial  accelerations  (Stokes’ 
slow-motion  assumption  consistent  with  linear  eddy 
dissipation  and  bottom  friction;  see,  e.g.,  Schlich- 
ting,  1968) 

(2)  Neglect  all  centrifugal  accelerations 

( 3 )  Neglect  vertical  Coriolis  force 

(4)  Neglect  vertical  dissipation 

(5)  Neglect  vertical  motion 

(b)  Single-layer  ocean  assumption  in  Equations  1  and  2:  i.e.. 

(1)  r  =  R  +  R,  but  dr  =  dz 

(2)  C(A ,*,r)«ff(X,0) 

C»(A,*,r)«ff(X,0) 

(3)  u(\,t;z,t)~u(K+,t) 

v(A,  0,  Z,  t)  ~  v(X,  4,  /) 
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It  may  be  pointed  out  that  all  assumptions  of  (a)  and  (b)  are 
well  justified  over  most  flow  areas.  In  particular,  the  strong  assump¬ 
tions  (b,  3)  are  realistic  because  the  tide-generating  potential  is  a 
body  force.  Moreover,  the  motion  is  fully  turbulent  and,  hence,  the 
averaged  velocity  profile  exhibits  only  a  very  thin  boundary  layer 
and  laminar  sublayer  (see,  e.g.,  Schlichting,  1968).  It  may  be 
mentioned  that  for  the  same  reason  the  condition  of  free-slip  at  a 
boundary  wall  (section,  “Lateral-Boundary,  Initial,  and  Final  Data”) 
seems  more  appropriate  than  the  no-slip  condition  used  in  laminar- 
flow  situations.  The  assumption  (b,  3)  is,  of  course,  much  less 
realistic  in  general  ocean  currents,  which  are  driven  by  surface 
pressure  and  wind  forces  and/or  by  density  variations  confined  to 
the  upper  ocean  layers. 

By  applying  assumption  (a)  to  Equation  lc,  one  finds  the  hydro¬ 
static  pressure 

P*  —  ~Gp, 
i.e., 

r 

p-p*=J  Gpdz  =  Gp(C  '  z), 

9 

or,  with  (section,  “The  Navier-Stokes  Equations  of  Averaged  Tur¬ 
bulent  Flow”  (a)-(g)) 

P*  =  P. 

r=;+*‘=c  +  r-r, 

(16) 

p=Gp(t  +  c-r-z)  +  p. 

With  Equations  7-10  and  1 6,  one  has 

=  G[v  ~  (**-!>')-*+ (r-f"  +  7*))  - Phr 


X 
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l.e.. 


q  ~  =  <7  (on?  -fit)-  P/p. 

P 


where 

<*  —  0.69, 
fi  =  0.90. 


(17a) 


(17b) 


In  the  following  equations,  it  is  convenient  to  introduce  the 
notations 


e+=h+/h,  e^h^/h. 


and 


H+  —  E+  —  Msin0/  ( 1  +  /icosd) 
so  that  (Equation  6a) 

A  -AH  ,  A  -AH  . 


(18) 

(19) 

(20) 

(21) 


Using  the  simplifications  (a)  and  (b)  and  Equations  17-21  in 
Equations  1  and  2,  one  arrives  at  the  reduced  equations  of  motion: 


Ut- 


R  cosd 


(«v  -  fit)y  +  20  v  sin  0  +  f\ 


v«  («>  -  0C)#  -  20  u  sin  d  +  ft. 


(22a) 


and 


U  +  (v  COS  0)#  +  U  COS  0  We  »  0. 


(23) 
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where  (Equations  3a,  b;  Sa-f;  and  21 ) 

“  R2  £s“0  t“xx  +  2IJ\U\  “«] 

+  ^2  iu<p<t>+  (H*  ~  tan0)u*  +  H+u  tan  0j 

+  Wcos*  ^g0  ~  2  tan  0)v\  ~  2^x  v  tan  0j,  (24a) 

>4 

f*  =  R2~cos~i  Kx  +  ^xVx  -  v]  +  ^v„ 

+  £2  [v^  +  (2#*  -  tan  0)V^  ] 

+  ^3^  K2«x+gx«)tan0  +  gx^J.  (24b) 

Under  the  single-layer  ocean  assumption  (b),  the  reduced  Equa¬ 
tions  22  and  23  may  be  integrated  over  the  instantaneous  ocean 
depth  (z"  —  zb  =  H  +  {)  while  observing  the  surface  and  bottom 
boundary  data  specified  in  the  section,  “The  Navier-Stokes  Equa¬ 
tions  of  Averaged  Turbulent  Flow.”  For  that  purpose,  it  is  useful 
to  introduce  the  “integrated”  velocity  components: 

U  (X,  0,  t)  =  J**  udz  ~uH —  (integrated)  east  velocity 

(25a) 

and 


V  (X,0,  t)  =  j*|J  ydz  ~  vff  =  (integrated)  north  velocity, 

(25b) 

where  the  term  “integrated”  may  be  omitted  when  no  confusion 
appears  possible. 

Using  the  simplifications  and  notations  made,  one  finds  the  help¬ 
ful  approximations  with  die  corresponding  replacements  (u-»v, 
U->V,X-*0): 


% 


Ocmh  TMms,  Pvt  L  Gfafevl  Occm  TMal  E»aHnt  177  f 

J*£  v«dz  ~  Ut  —  v'z'r  +  kVi  **  (/»,  (26a) 

JV  uxdz  =  UX-  Vz*x+  *VX~  Ux  ~  tix  U.  (26b) 

and 

J?  «xx* ~UXX-HXUX  +  (J?i -  J7XX)(/  - »x«l 

(26c) 

The  bottom-boundary  conditions  imposed  in  the  section,  “The 
Navier-Stokes  Equations  of  Averaged  Turbulent  Flow,”  (j)  and 
(k)  assume  the  following  approximate  forms: 

^=S^?“*+TV+W‘-  <27> 

— T*>x  =-—  U  •“ A  ( — +^^  +  11*^  (28a) 

p  HU  A\R*co^^  R*  (Z8a) 

and 

7*=i^(*wv+i**+4  <28b> 

With  Equations  25-28  and  die  surface-boundary  conditions 
in  die  section.  The  Navier-Stokes  Equations  of  Averaged  Turbu¬ 
lent  Flow”  (h)  and  (i),  the  proposed  integration  of  Equations  22 
and  23  is  easily  curled  out  and  yields  the  following  “continuous 
ocean  tidal  equations”  (COTEs  with  e  *  »/2  -  ♦  =  colatitude): 


+  t/„  +  <cot6  -  V,  +  H,  )l/,j  (29.) 

+  F  ~  (2  cot  e  +  i?# )  Kx  J  +  2nF  cos  e, 


t 
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+  +  +<“‘«-ff.+2H.)'/,j  (2»b) 

+  jftffr-Q  [ 2 cot 6  t/x  -nxUt  -  (cot  $  -  H,)HXU)  ~  2nVco&6, 

and 


Jlsme&  +  Vx  -  (Krin  0)#  *0,  (30) 


where  «  =  0.69,  p  =  0.9,  and 


ffx  _  ft*  _  ^  (Hr#  -  U9  )  +  (H,  +  Hf )  cot  6, 

(31a) 

H* s  +  *w  +*#  («*9 -  »♦  +  2ff# ),  (31b) 

ff,  =  (6 -*  A.)  \ 

J  (32) 

B0  =*I7#  +  *coe0/(l  +pm  0),  ] 

and 

,4  »|LH(A,e)0  +  *rine),B*  b^sinS.  (33) 

It  may  be  mentioned  that  for  «  =  0  =  1  and  A  =  B  =  0  the  ocean 
tide  equations  (Equations  29  and  30)  reduce  to  die  considerably 
simpler  classical  Laplace  tidal  equations.  Evidently,  the  complete 
COTEs  require  second  derivatives  (Hu.  H„)  of  the  bottom  topog¬ 
raphy,  which  can  be  assumed  to  exist  without  placing  major  re¬ 
strictions  on  the  realistic  features  of  the  bathymetry.  It  is  interesting 
to  note  that  those  second  derivatives  in  H*  and  H*  act  as  modifying 
bottom-friction  terms.  A  ridgetike  ocean  floor  (say,  Hu  >  0) 


always  adds  to  the  bottom-friction  terms  (U,  V)  •  B/H,  while  a 
valley  (H  u  <  0)  always  diminishes  bottom  friction. 

Second-Order  Arctic  Tides 

The  COTEs  (Equations  29  and  30)  become  singular  at  the  North 
Pole.  For  the  intended  numerical  procedure,  it  is  therefore  advan¬ 
tageous  to  seek  an  approximate  analytic  solution  around  this  sin¬ 
gularity  that  can  be  matched  together  with  the  numerical  solution 
at  some  appropriate  colatitude.  In  fact,  a  unique  “second-order 
arctic  tide”  solution  can  be  determined  for  all  throe  species  of  tide- 
generating  potentials  listed  as  Equations  12,  13,  and  14,  provided 
the  ocean  depth  around  the  North  Pole  is  assumed  constant. 

For  constant  depth,  H  =  Ho,  constant  eddy  viscosity,  A  =  Ao, 
and  constant  bottom-friction  coefficient,  B  =  Be,  the  COTEs 
(Equations  29  and  30)  assume  die  simpler  form: 


+  Ratin’©  £/  +  sin0  (sin  0  U$  )t  -2cos0Kx], 

(34a) 

GHo  Bo 

=  W-*v),  -2flI/cose~grK 


and 


+ 


Aq 

R*  sin’e 


-  K  +  sin  0  (sin  0  Vt  )$  +  2 cos  0  £/x]. 


(34b) 


fc  +  j^Q[tfx-(Psine)  J=0.  (35) 

The  forcing  equilibrium  tides,  n  (Equations  12,  13,  and  14), 
may  be  written  in  the  unified  complex  form 


v  *?  (A,  0)  e*. 


(36) 


where 
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9)  = 

With  the  substitution 


ffsin*0e*\  (*-  2)  (37a) 

sin  26  e*S  (»  =  1 )  (37b) 

K( 3  sin3  9  -  2)  (37c) 


<v,  C,  U,  V )  «  (7>  f, Z7„  iT)eV  (38) 


one  arrives  at  the  following  three  complex  differential  equations  in 

(A,  ©): 

+2oT'cos6+i^1' 

(^-P  +  metsinei;,  >,  -2/cmSVJ, 

(39a) 


~V=^j~-(<*r~0Z)9  +2aUcosQ  +  ijj0V 


iAo 


-  fF*  -V~+  sin0  (sin 6  F#  ),  -  2/cos9l7x], 

(39b) 


and 

•f=«4T|iF‘  +  (p,“!)>1' 

(40) 

The  reduced  Equations  39  and  40  may  be  solved  for  regular 
solutions  by  power  series  in  sin  6  (essentially  polar  distance)  of  die 
form: 


For  *  ~  2  and  »  =  0, 


7®  fr(A)  +  *i(a)  sin 6  +  6(A)  tin*  6  +  . . .. 
V  ■  tf.(A)  +  £/,<A)  till  e  *  Ut( A)  tin*  e  +  . 


*  •  » 


and 

F  *  cot  6  [Fo(A)  +  ^i(a)  sin  6  +  F*(a)  sin*  6  + 


and  for  *  =  1, 


Z=  cos  8  U*(a)  +  ti(A)  sin  0  +  6(A)  sin*  0  +  . . .  ],  \ 
ff  =  cotfl  [£/t(A)  +  f/i(A)  sin  6  +  I/,(x)  sin*6  +  . .  .J 


and 


X42) 


F=F*(a)  +  F,(a)  sind  +  16(A)  sin2  0+ _  ) 


By  truncating  the  series  expansions  (Equations  41  and  42)  after 
die  second  order  in  sin  0  and  substituting  these  truncations  into  the 
Equations  39  and  40  up  to  the  same  quadratic  power,  one  arrives 
after  some  lengthy  but  simple  algebra  at  the  following  unique 
“second-order  arctic  tidal  approximations.”  For  semidiurnal  tides 
(v=s  2), 


1  —  J?sin*6  ei<ix+*l\ 

C  =  ti^sm*©  e<,*x+»t', 

U  —  —2K.  <r  it  tin  0  e‘,*x+»l), 

and 


M43) 


F  =  -flr®rttm*e**,3*+'\ 


where 


K  =  «GffO:/[6/SG«o  +  o**(  2n  -  «r)  +  fr{6A*  +  B*R*/Ho) ) ; 

(43a) 

and  for  diurnal  tides  (»s  1), 
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r,  =  K  sin  26  1 

t  =  3^i  sin  2fi  l  (44) 

U  (&>  +  6^3  sin2  6)  cos  9  e*'*+*“, 

and 

V  =  i  [TCa  +  2(K*  +  Jti  *R)  sin2  9]  «*'fc+*\ 

where 


^  KhKMKs*- Ku) 

Kt~Ka( KtsKsi  -  KnKn)  -  Ku(KnK* *  -  Jhb)' 

=  -  Xa(*is  -  JTu  JCm/ITm)], 

A.X2 

Xx  —  “RAlm/JCji; 


(44a) 


and 

JTu  =  (60GHo  -  «*/?*)  +  /«r(64.  +  BaRVffo), 

ATxa  =  -flR, 

tf.i »  K(6n  -  <r)  +^~(  l2Ao  t  B*R*/Ho), 

Ki\  ~  6jJGH o  +  4ivAo, 

Kaa  =  /*(*  -  20)  -  £  (2WU  +  BoR*/Ho)f 

Kit  -  16iAa/R, 

Kti  *  2odR\ 

Ku  *  R(2Q  -  3<r)  +  3 ~ ( 12^o  +  BoR'/Ho), 

Kn  =  2*GHoK. 

For  long  period  tides  (v  =  0), 

7  =  -KJ2  -  3  sin^)***, 

{=  -2Ki(2  -  3  sin20)e^* 
i/  =  A?j  sin  9  «***, 

and 


>  (44b) 


\ 


j  (45) 


K  =  -  i%i<rR  sin  20  «*•', 
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where 

K,  =  4«OffV[«rK*  ~  i(2Ao  +  BoRVHo)  ], 
K,  =  ~KiK*, 


and 


(45a) 


Ki  =  -3aGHoK  [  (6/3GHo  ~  o*R*  +  ORK*) 

+  i<r(6A*  +  B*R3/Ho)]. 

It  is  important  to  note  that  the  surprising  uniqueness  is  achieved 
by  requiring  a  regular  integral  at  the  North  Pole  (Equations  41  or 
42)  and  by  truncating  the  series  of  the  solution  and  the  differential 
equations  after  the  second  power  of  sin  0.  Of  course,  without  the 
second-order  truncation,  uniqueness  can  no  longer  exist  Undeter¬ 
mined  coefficients  become  available  to  satisfy  prescribed  boundary 
data,  say,  at  distant  continental  shorelines. 

In  the  present  global  tide  model,  the  arctic  tides  (Equations  43, 
44,  and  45)  will  be  considered  valid  up  to  1°  colatitude  (the  first 
land  occurs  at  colatitude  7°;  see  Bathymetric  Tables  in  Schwiderski, 
1978a).  For  colatitudes  2°  and  3°,  a  linear  interpolation  will  be 
used  to  match  the  polar  solution  with  the  numerical  solution  com¬ 
puted  south  of  3°  colatitude. 

It  is  interesting  to  observe  that  if  the  Coriolis  force  and  eddy 
dissipation  are  neglected  (ft  =  0,  Ao  =  0),  then  the  second-order 
arctic  tides  (Equations  43,  44,  and  45)  become  exact  global  tides 
with  the  constants 

=  aGHoK/[  (60GHo  -  o*R2)  +  iaBJP/H*], 

—  2K,  Ts  —  ~2oRK,  —  0, 

kr  =  -3X,  K*  =  0,kt  =  0. 

From  Equations  43,  44,  and  45,  one  concludes  that  all  second- 
order  arctic  tides,  {,  vanish  at  the  North  Pole  with  the  same  order 
as  their  corresponding  driving  equilibrium  tides  7.  Hence,  only 
long-period  tides  exist  at  the  North  Pole. 
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Derivation  of  Discrete  Ocean  Tide  Equations 
The  l9byl°  Graded  Grid  System  and  Bathymetry 
With  the  exception  of  Antarctica  (south  of  colatitude  6  =  168°), 
the  entire  (ocean  and  land)  area  of  the  globe  is  covered  by  a  1°  by 
1°  grid  system  that  is  “graded”  toward  the  poles.  Each  spherically 
rectangular  mesh  cell  Sm.n  is  bounded  on  the  east  and  west  by  longi¬ 
tudes  An,  =  m°  andAm_*  =  (m  -  m)°,  respectively,  and  to  the  south 
and  north  by  colatitudes  9n  =  n°  and  9„-i  =  (n  -  IV,  respectively, 
so  that 


Sm.n  —  ( ff! 


mV 


(n-iy 


where 


m  ~  ii,  2m»  • . .  ,(360-»0) 

and 

n  =  1,2,. .  .,168. 


(46) 


(47) 


The  “grading”  parameter  m  =  is  defined  by: 

m  =  1  for  n  =  30  to  150 
**  =  2  for  n  =  15  to  29  and  n  =  151  to  168 
/*  =  4  for  n  —  8  to  14 
a»  =  8  for  n  =  1  to  7 


(48) 


The  grading  of  the  network  toward  the  poles  js  necessary  in  order 
to  maintain  a  more  uniform  mesh  area  for  higher  accuracy  and 
stability  (section,  “Stability  Analysis”)  of  the  discrete  tide  model. 
In  fact,  the  grading  equations  (Equations  48)  have  been  chosen  in 
such  a  way  that 


Main  n°  >  1/2  for  n  =  4.  8,  15,  30, 150, 165;  (49) 


Le.,  the  southern  mesh  size  remains  larger  than  half  the  equator 
mesh  size.  This  important  condition  is  violated  for  n  =  1°,  2°,  3° 
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and  n  =  166°,  167°,  168°.  However,  in  the  preceding  section  it 
was  pointed  out  that  the  numerically  discrete  tide  model  begins  at 
colatitude  9  =  4°.  For  colatitudes  9  =  3°  and  2°,  the  numerical  so¬ 
lution  will  be  matched  by  linear  interpolation  to  the  second-order 
arctic-tide  approximations  derived  in  the  preceding  section,  which 
are  assumed  to  be  useful  up  to  colatitude  9  =  1°.  As  will  be  shown 
in  the  section,  “Stability  Analysis,"  the  slight  violation  of  the  con¬ 
dition  in  Equation  49  at  the  three  southern  colatitudes  is  not  so 
severe  as  to  affect  the  stability  characteristics  of  the  model.  In  any 
case,  no  need  was  apparent  to  justify  an  additional  grading  step, 
which  is  accompanied  by  an  unnecessary,  extra  computational  ef¬ 
fort  (next  section) . 

In  the  global  network  defined  above,  land  and  ocean  mesh  cells 
are  distinguished,  respectively,  by  zero  or  nonzero  depth  data  H.,, 
which  will  be  assigned  to  each  cell  below.  The  “mathematical” 
boundaries  of  the  oceans  follow  in  an  obviously  zigzagging  fashion 
the  mesh  lines  of  boundary  cells. 

The  ocean-depth  data  collected  by  Smith  et  al.  (1966)  were 
rearranged  and  linearly  interpolated  to  fit  the  new  1°  by  1°  graded 
grid  system  described  above.  The  original  data  bank  had  to  be 
corrected  for  obvious  errors  in  continental  and  oceanic  labeling, 
ocean  and  land  signs,  shorelines,  and  some  exponents.  All  land 
elevations  were  set  to  zero,  including  all  depth  data  less  than  50  m. 
Furthermore,  the  following  meshwise-disconnected  border  seas 
were  excluded  from  consideration  by  assigning  zero-depth  values  to 
their  corresponding  mesh  cells:  the  Baltic,  Kattegat,  Irish,  Mediter¬ 
ranean,  Red,  Japan,  Sulu,  and  Ceram  seas;  the  Hudson  and  Korean 
bays;  and  the  Chihlian,  Persian,  and  Californian  gulfs. 

All  depth  data  Hm,»  are  considered  representative  for  the  center 
of  the  cell  Sm.D;  i.e.,  for  m  =  n,  2/x. .  .,360  and  n  =  1,2,. .  .,168, 

H».b  =  H(a«,9.)- 

where  (m  —  #i/2)#  and  0»  =  (n  —  Vi  )*\ 

The  assigned  minimum  and  maximum  depth  values  were,  re¬ 
spectively, 
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The  averaged  North  Pole  deoth  was  found  to  be 
Ho  =  3,600  m. 


(51) 

(52) 

(53) 


which  is  needed  to  compute  second-order  arctic  tides  (section, 
“Second-Order  Arctic  Tides”)  for  colatitude  line  n  =  1. 


The  Discrete  Ocean  Tide  Equations  (DOTEs) 

Following  essentially  the  hydrodynamical  numerical  method  of 
Hansen  ( 1966)  and  Zahel  (1970),  the  COTEs  (Equations  29  and 
30)  may  now  be  converted  to  an  explicit  finite-difference  analog 
called  “discrete  ocean  tide  equations”  (DOTEs).  In  agreement 
with  the  graded  grid  system  defined  in  the  preceding  section,  it  is 
convenient  to  introduce  the  following  notations: 

A0  =  ir/36O=  V4°  mesh  size,  (54) 

AA  =  (jAQ, 

L  —  2RA0, 


At  =  time  step  to  be  specified. 


(55) 


In  a  first  step,  all  spacial  derivatives  of  Equations  29  and  30  are 
replaced  by  divided  central  finite  differences,  using  a  Richardson 
(1922)  staggered  scheme.  Accordingly,  for  a  fixed  oceanic  mesh 
cell  Sn.n  (Equation  46;  j  =  1,  2, 3 . . . ) 


U-.i  =U(AB(u),9„(U),jAt), 

vL+.i  =v(AB1(v),e„(v),jAt), 

=<(A.(0,9.U>,  j*t) 


(56) 


are  computed,  respectively,  at  the  u-point  of  S».« 
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x»(m)  =  2(m  -  m)ax,  6»(u)  —  (2n  -  1)a6,  (57a) 

v-point  of  Sb.b 

a*(v)  =  (2m  -  S»(v)  =  2n&6,  (57b) 

and  {-point  of  Sm.D 

A.(()=Mv),6.  U)=e.(«),  (57c) 


Undefined  data  are  specified  by  arithmetic  means  of  defined  values. 
For  example,  needed  U-data  at  v-points  are  defined  by  averages  of 
U-values  at  adjacent  u-points. 

In  the  second  part  of  the  differencing  process  the  resulting  equa¬ 
tions  are  integrated  in  time  over  a  single  time  step,  At,  by  using  the 
following  average  integration  rule  (U-»V-»{) 

f  U(t)dt  =  At  [*Ui+1+  (l-*)U')t  (58) 

j  if 


where 


U*  ~  U (tj) ,  tf~  (/  1 )  At,  /  =  1,2 . 

and  «  =  some  differencing  parameters,  which  usually  satisfy  the 
restriction 


(59) 

As  will  be  seen,  the  resulting  discrete  ocean  tide  equations  (DOTEs) 
are  very  sensitive  to  the  choice  of  the  differencing  parameters,  «, 
and  the  time  step,  At.  In  fact,  depending  on  the  chosen  values  of  * 
and  At,  the  DOTEs  may  be  stable  or  unstable  (section,  “Stability 
Analysis”)  for  any  specified  values  of  eddy  viscosity.  A,  and  bot¬ 
tom-friction  coefficient,  B.  Moreover,  different  values  of  «  for 
Equations  29  and  30  and  for  die  various  point  values  of  U,  V,  and 
(  can  be  chosen  so  that  the  resulting  DOTEs  may  become  explicit 
or  implicit.  Considering  the  complexity  of  the  ocean  basin  and  the 
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large  number  of  oceanic  mesh  cells,  an  explicit  form  for  the  TXDTEs 
is  here  preferred. 

After  carrying  out  the  described  integration  with  the  three  (ob¬ 
vious)  differencing  parameters  *  =0,  *  and,  5,  one  arrives  at  the 
following  explicit  DOTEs.  (For  a  more  detailed  derivation  see 
Schwiderski,  1978b.) 

[l  +  «4S,]  =  A„:.Hsin^I(2y-  1) 

-  A*  „  COS  (2/  -  1)  +  Am,n  [&-*.«  ~  &*] 

+  [1  -  (1  -  K)A,1„]  +  Afn.nUjm^n 

A  Am,nUin-ti.n  A  A]aM  t/m.n  +  l  +  Am,n^m.n- 1 

+  A®,„  [(Vi,.„  +  Vi,,.,)  -  (Wm.M(B  +  Vi,.^.,)] 

+  A  J®.„  [(VJm.„  +  Vi,.*-,)  +  (F4_m.«  a  (60a) 

[i  +• « c,  ]  v,:v  =  cos  ^  <2 j  - 1) 

+  ««.  sin  (2j  -  1) 

+  [U..I  -  a„]  +  [l  -  (1  -  *>«.*,,,]  Vl„ 

A  A  V/n-/ii,n  J  A  5,®lH  Vm.n  +  l 

+  R 7  VJ 

'  **«  ,rt  r  M  ,rt  -  I 

A  B,ln  [C/,i+M>n+I  -  £/>,„]  A  Bm9,„  [f/,i+M.n  - 
+  B"n  [Vi  +  f/i,,  +  Ui  +M.H  A  £4.,+.] ,  (60b) 


and 


V* 1 


c,, + « K'  +  WJ  - 

+  (I  -  [c.'d/i,  -  e/i — ) 

+  c;v^-c;v>].  (W) 


% 


4 

i 
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The  coefficients  of  the  DOTEs  (Equations  60  and  61 )  are 

AmJt  =  Am,n  cos  2  v{m  -  PL)  Ad, 

A£,„  =  Am,n sin  2 v(n  -  pl)  Ad, 

Am.n  =  f3AtGH{u)Yn(u)IL, 

m.n 

Am.n  =  2a  4^  *MW)  [<»„(«)//(«)  -  tf(i/)] 

^  m.n  m.n 

+  bAtlYn{u)H(u), 

m.n 

A'm.n  =  fl  4 ”  »^(")r«  {u)H{u)  [l  +  //(«)] , 

A*m.„  =  a  ±t-Mu)Y,?(u)H(u)  [l  -  //(«)], 


m.n 


m,n 


Al.n  =  b~  H(U)  [*„(*)  +  ^  (2  +  r„(«)) 

m  ,n  ^ 

*cos(2/i  -  1)  Afl] , 

4L.  “  *>  •¥  H(u)  [«M«>  -  ^  (2  +  r,(«)) 

^  mm  ^ 

•tos(2n  -  1)  A0] , 

K,  =  O  4 -  I »w<«)  [*.(«)#<■<)  -  iiM  (3  +  2T,(»)) 

m.n  m.n  * 

. cos (2 n  -  1)  A0] , 


and 


(62a) 


A"‘«  =  -  «  4 Lr»WH(u)H(u)  [*Au)ft(U)  -  UM 
cos(2«  -  l)A0]  +^(lCos(2n  -  l)  SO, 

where 


A 


nt  Jt 


~  4  ~j?KH(u)  sin  ^ 

w.i>  2 


sin  (2n  -  1)  AO, 
cos  (2/i  -  l)A0, 


0 


and 

Liu)  =  1,V  sin  S,(w);  ^„(«)  =  ^-(1  +  t/TM(t,)), 
<*^(«>  -  i  +  r;»; 


B,l,.n  ~  B,,,  ,,  cos  1/(2 m  -  ji)  £0, 


#«.«  -  -  Bm„  sin  v(2m  ~  fi)  AO, 

B,X,,  -  &AtGH{\)iL, 

m,n 

B>ln  ®  2"  X  <MV)[<4,(V)//(V)  ~/7(V)] 

^  m,n  fn.H 


i-  b  AtlY„{\)H(\), 

mjt 


B*  «  '■ 

•*4i»  m 


“  «  ■—  **W(v)tf(v), 

m.w 


K.„  *  «  ~  *<v>  K(V)  ( 1  -  ^(v) ) 

*-*  Hi  .H  mm  mm 


\ 


> 


>  (62a) 


v  *  2 

*>  *  1 
»/  =  0 

(63a) 


(62b) 
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+  Ji|i(3  +  2C(v))  cos  2/7  Ad] , 
K.„  =  a  ~  H{\)  [«Mv)  (1  +  tf(v)) 

**  mj9  m,n 

-  (3  +  2r„(v))  cos  2nA0], 

Sm.H  ~  «  4^  ^(V)0„(V)//(V) 

m,« 

[2#iA^r„(v)  cos  2nA6  -  tf(v>] , 

m.n 

K.n  =  «  —  r„(v)i|ifl(v)//(v) 

*-»  m.n 

[2/LiA0r„(v)  cos  2/tA0  +  //(v)] , 

m.n 


/ (62b) 


#m°«  *  ~  a  ^  r„(vW„(v)H(v)H(v) 

*-  mm  m.n 

[/xAfll^v)  cos  2nAB  +  tf(v)]  -  cos  2nAB, 

m,n  £ 


where 


2  sin  An  AS,  v  =2 
cos  An  AS,*  ,  >>  =  I 
6  sin  4/t  A0,  k  =  0 


and  >' 

j 


*  *  4ff  KH^]%in£r~ 


192 


Erast  W.  Schwiderski 


r.(v).  =  l/fisin  0,(v);  0„(v)  =  -l(l  +  1/C(v)). 


(63b) 


< w.,(v )  *  1  +  T.2(v); 

c««-  —-Guo, 


q*>  *  M  ^  f»  sin  2nA0. 


>  (64) 


and 


^  ..  A/ 


C®  =  m  r„(«)  sin  2(n  -  1)A0. 

In  these  coefficients,  the  following  depth  functions  were  used  at 
u-points: 

H(u)  = 

/n./t 


"<“> -  iWm <*> -  h <">  ]• 

m,n  •filyU)  m-n,n  m  m.« 


MI.II 


#(«)  38  -TJTTut  [H  (")  ~  «  («)  ]. 

*»*  mji-i  «.*+l 


#(«)  - 


//(«) 


<•»*(«)  -  2C2(W)  [H\u)  +  (M*)2] 

DM 

+  cos  (2n  -  X)A0  [(3  +  21 *(«))#(«) 

-  ^A0r,(i/)  cos  (2 n  —  1)A0] 

-  •i[r,»(«)  (h  («)  +  //<«)  ) 

d*  m  +  M*  m~u*H 


+  («(«)  +  « («>  )]; 

ntji-t-1  m,n-l 


(65a) 
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and  at  v-points 

H(\)  «  i  [HmJ,  +  HM'H+l], 

MJi 

»(*)  =  4HM  t"<V)  ”  W(V)  1 

w »n  H/l  w+M 


m,n 


W<v)  "  3777^  IH<v)  "  «(v)  ] 

m.n  4/1  (V)  m.n  —  1  m.n  +  1 


£(V)  =  //(V) 


m,w 


IN, ft 


o»„(v)  -  2[tf2(v)  +  (/*A0)*i;*(v)] 

ftl.ll 

+  /i.A#ih,(v)/^(v)  [3  +  £,(v)]cos  2nA0 
-irn*(v)(//(v)  +  #(v)  ) 

^  ftt  +  M*N  IN—  JMt 

+  («(v)  +  «<v)  )]. 

m,n+ 1  m.n-1 


(65b) 


Obviously,  due  to  die  introduction  of  the  novel  depth-dependent 
eddy  viscosity  (Equation  6)  the  COTEs  and  the  DOTEs  (Equa¬ 
tions  29,  30  and  60,  61)  are  considerably  more  involved  than  die 
corresponding  equations  used  by  Zahel  (1970,  1973,  1975)  and 
Estes  (1975).  This  results  in  drastically  increased  computer  time, 
memory,  and  cost. 

For  *  =  Oand  «  =  1,  the  finite-difference  scheme  coincides  with 
the  technique  used,  e.g.,  by  Hansen  (1966),  Zahel  (1970,  1973, 
1975),  and  Estes  (1975,  1977).  Extensive  exploratory  computa¬ 
tions  were  carried  out  by  the  author  with  numerous  *  and  «  values 
within  the  ranges  0  1.5  and  0.5  <  *  <  1.  The  computations 

produced  no  drastic  differences,  provided  the  eddy  and  bottom- 
friction  coefficients  a  and  b  and  the  time  step  at  were  suitably 
chosen  within  their  respective  stability  constraints  (section,  "Stabil¬ 
ity  Analysis") .  Finally,  it  was  decided  to  use  the  values 
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K  —  K  —  1 


(66) 


because  they  yielded  a  preferable  stability  and  seemingly  best  re¬ 
sults.  Moreover,  this  deviation  from  the  Hansen-Zahel  method 
became  most  significant  for  the  hydrodynamical  interpolation  of 
empirical  tidal  data  described  in  Part  II  (Schwiderski,  1979b). 
Indeed,  this  novel  technique  uses  the  special  property  that  for  *  =  1 
the  bottom-friction  coefficient,  b,  which  enters  only  in  A4  and  B4 
(Equations  62)  becomes  essentially  a  scaling  multiplier  of  Um  and 
in  Equations  60a  and  60b.  Thus,  together  with  t  =  1  in  Equa¬ 
tion  61,  the  bottom-friction  coefficient  can  easily  be  adjusted  locally 
to  match  more  closely  prescribed  tidal  data. 

With  *  and  *  specified  by  Equation  66,  the  DOTEs  (Equations 
60  and  61 )  still  contain  the  parameters  a,  b,  and  At,  which  remain 
at  one’s  disposal  within  their  respective  stability  ranges  (section, 
"Stability  Analysis”).  They  will  be  utilized  to  achieve  best  results 
by  trial-and-error  computations.  The  DOTEs  (Equations  60  and 
61)  can  be  applied  to  all  oceanic  mesh  cells  Sm  D  with  m  =  m,  2m. 

. . .  ,360  and  n  =  4,5 . 168  sweeping  across  the  globe  from  n  = 

4  to  n  —  168.  This  procedure  can  be  executed,  provided  suitable 
initial  and  lateral  boundary  data  (next  section)  are  prescribed.  At 
colatitude  line  n  =  4,  the  numerical  solution  is  matched  to  the 
second-order  arctic  solution  (section,  “Second -Order  Arctic  Tides”) 
by  the  linear  interpolation  (m  =  m,  2m,  . . .  ,360) . 


}  (U-+V-+ C) 


(67) 


where  U£*  V£|  and  are  computed  by  Equations  43,  44,  or 
45.  For  colatitude  lines  n  =  7,  14,  29,  150,  specially  corresponding 
data  (U,  V,  ()  on  n  =  8,  15,  30,  151  (see  Equation  49)  are  de¬ 
fined  by  linear  interpolation.  Vice  versa,  for  n  =  8,  15,  30,  151, 
specially  corresponding  values  (U,  V,  {)  on  n  =  7, 14, 29,  150  are 
also  defined  by  linear  interpolation. 
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Lateral-Boundary ,  Initial,  and  Final  Data 
In  order  to  complete  the  ocean-tide  model,  the  DOTEs  (Equations 
60  and  61 )  must  be  supplemented  by  suitable  lateral-boundary  and 
initial  values.  In  turbulent  flow  situations,  the  mathematical  boun¬ 
dary  conditions  usually  preferred  are  (a)  no-flow  across  the  ocean 
shorelines  and  (b)  free-slip  along  the  ocean  shorelines.  It  is  clearly 
at  this  point  that  the  great  attractiveness  of  Richardson’s  (1922) 
staggered,  finite-difference  scheme  (preceding  section)  manifests 
itself  in  the  practical  simplicity  with  which  the  no-flow  and  free- 
slip  (or  no-slip)  boundary  conditions  can  be  worked  into  the 
model.  In  fact,  if  &■,.  is  an  oceanic  boundary  cell,  then,  by  defini¬ 
tion  (section,  “The  l°xi°  Graded  Grid  System  and  Bathymetry”), 
the  zig-zagging  mathematical  boundary  lines  follow  only  mesh  lines 
and  pass  only  through  velocity  (u  and/or  v)  points.  The  no-flow 
condition  (a)  is  implemented  by  declaring  at  those  points 

l/^lm=  0  and/or  V’*\  =  0,  (68) 

respectively.  If,  say,  the  v-point  of  S...  is  oceanic  and  the  v-points 
of  S  and/or  S  are  terrestrial,  then  the  free-slip  condition 
(b)  is  satisfied  by  reflectively  setting 

Vi+1  =+V>+1  and/or  F'+I  =  +V*'\  (69) 

respectively.  (If  the  no-slip  condition  is  imposed,  then  the  (+) 
signs  in  Equation  69  must  be  changed  to  (  - )  signs. ) 

The  mathematical  boundary  conditions  (a)  and  (b)  were  ap¬ 
plied  by  Hansen  (1966),  Zahel  (1970,  1973,  1975),  Estes  (1975, 
1977),  and  others.  In  the  present  model,  both  the  free-slip  and  no¬ 
slip  conditions  were  tested.  The  computer  experiments  indicated 
slightly  better  results  for  the  free-slip  condition,  which  was  then 
adopted.  Physically,  this  condition  is  more  plausible  in  turbulent 
flows  which  have  only  very  thin  boundary  layers  (see,  e.g.,  Schlich- 
ting,  1968).  Furthermore,  the  free-slip  condition  is  consistent  with 
die  bottom  boundary  assumption  specified  in  the  section,  ’The 
Continuous  Ocean  Tide  Equations  (COTEs)”  (b,3). 
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The  construction  of  the  Ms  tide  was  started  at  j  =  1  (ti  =  0)  with 
an  ocean  completely  at  rest;  i.e.,  with  the  initial  values  (m  =  m. 
2/*,.  .  .360;n  =  l,2,...168) 

U1  =  V1  ={»  =0.  (70) 

m.n  an  m.n  w  / 

The  computations  were  carried  over  a  prescribed  number  of  time 
steps,  j  =  J  (mostly  a  quarter  period),  and  then  printed  for  inspec¬ 
tion  of  the  results.  With  or  without  program  and/or  parameter 
changes,  the  computations  were  restarted,  using  the  latest  or  any 
earlier  taped  output  instead  of  the  initial  values  (Equation  70). 
Occasionally,  it  was  beneficial  to  speed  up  an  unwanted  slow  decay 
of  transient  eigenmodes  by  “negatively”  averaging  the  output  data 
of  half  a  period  time  difference;  i.e.,  by  setting 

UJ  =  Vi  UJ  -  UJ-'  ,U-+V->(,  ill) 


where  <rAtJ  =  w.  This  simple  procedure  diminishes  all  undesirable 
eigenmodes  of  lower  frequency  than  the  forced  frequency,  <r,  and, 
similarly,  also  most  higher  frequency  modes.  The  negatively  aver- 
aged  data  (Equation  71)  represent  obviously  improved  initial  data. 

Two  output  tidal  elevations  2J-time  steps  apart  (mostly  a  quarter- 
period  apart), 


and 


=  *«,.  cos  [•*(/  -  27) V.],  j 


were  used  to  compute  the  tidal  “amplitudes” 


u*  =  [<\ + ni* 


(72) 


(73) 


and  “phases” 


!■,.»  =  arc  tan  y/x  (=  Ofer*  =  y  =0), 


(74) 


% 
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where 

it  =  (&?+  CU )/  2  cos  oA/7, 

(75) 

&=(t^-fiM.)/2«noAi7; 
y  =  ii  sin  <rt4(J  -  7)  -  is  cos  »A/(/  -7), 

and 

jr  =  fi  cos  <rA/(/  -  7)  +  6  sin  ®A/(/  -7). 


In  the  finished  product,  i.e.,  when  all  (unforced)  transient  eigen- 
modes  have  satisfactorily  decayed,  die  amplitudes  (Equation  73) 
and  phases  (Equation  74)  become  essentially  independent  of  time 
and  undergo  no  further  variation  of  significance  after  continued 
computations. 

In  the  MHide  model,  the  computed  convergence  toward  the 
steady  state  of  the  amplitudes  and  phases  was  found  to  be  gener¬ 
ally  oscillating.  So  the  integration  process  could  safely  be  termi¬ 
nated  when  the  amplitudes  and  phases  over  most  ocean  areas 
varied  by  less  than  2  cm  and  3°,  respectively.  The  convergence  was 
slightly  less  complete  in  coastal  cells  where  the  tidal  elevations  are 
extremely  large  and  vary  rapidly  from  degree  to  degree. 

In  order  to  follow  the  convergence  of  the  computer  program 
more  closely,  the  squared  tidal-amplitude  sum 

360 

*1=  1  .  ....  (77) 

w-m  . 

was  computed  and  printed  for  each  fixed  colatitude  line  n  = 
4,5, . . .  ,168.  To  compute  the  amplitudes  in  Equation  77,  Equa¬ 
tions  73  and  75  were  used  for  27  -  1  and  J  =  j  +  1 ;  i.e.,  for  two 
consecutive  time  steps  in  connection  with  Equation  61.  In  this 
measure  (Equation  77),  die  convergence  was  carried  to  almost 
three  significant  figures  for  all  n. 
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Stability  Analysis 

A  rigorous  stability  analysis  of  the  homogeneous  DOTEs  (Equa¬ 
tions  60  and  61)  is,  of  course,  not  possible.  However,  under  the 
assumption  of  constant  coefficients  A,  B,  and  C,  the  simplified 
DOTEs  possess  Fourier- type  eigensolutions  (Equation  88)  that 
permit  a  local  stability  analysis  of  the  difference  system.  As  is  well 
known  (see,  e.g.,  Richtmyer,  1957),  such  a  local  stability  analysis 
produces  stability  limits  that  are  usually  sufficient  for  computational 
purposes.  Indeed,  computer  experiments  showed  that  the  stability 
limits  so  derived  below  were  scrupulously  binding  for  the  success 
of  the  integration.  The  following  analysis  is  an  expanded  version 
of  die  investigation  presented  by  Zahel  ( 1970) . 

In  detail,  the  following  simplifications  may  be  assumed: 

(a)  b  -  0;  i.e.,  no  bottom  friction 

(b)  o  =  0;  i.e.,  no  Coriolis  force 

(c)  For  an  arbitrary  but  fixed  mesh  cell  (see  Equations 
62-65), 

T  =  r,(u)  =  r„(v)  =  1/psin  8  =  locally  constant, 

*  =  *„00=Mv)  =  **(l  +  l/T), 

<0  =  0>«(u)  =  <*>»(v)  =  1  +  r*; 
and 

H  —  H(u )  =  H(v)  =  locally  constant, 

7i(u)  =  H(u)  =  H(u )  =0,  (u->v). 

m, m 

It  may  be  mentioned  that  the  assumptions  (a)  and  (b)  have  been 
made  in  order  to  display  more  clearly  the  most  important  stability 
characteristics  of  the  DOTEs  that  are  due  to  eddy  dissipation  and 
to  the  differencing  parameters,  *  and  It  is  relatively  easy  to  show 
that  the  bottom  friction  is  always  stabilizing,  while  the  Coriolis 
force  (in  die  present  differencing  scheme)  is  slightly  destabilizing. 

For  the  following  derivations,  it  is  helpful  to  introduce  some 
specified  reference  values  9r,  Hr,  (consistent  with  Equations  78  and 
79)  r„  <(>t,  <*,  and 


(78) 

(79) 


t 
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Finally,  the  following  relative  quantities  may  be  introduced: 


1 1/2 
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(80) 


J  —  <p/<pr,  ft)  =  ft)/ (Or,  h  =  H/Hr 
r  —  A// A/r,  (  =  a,' Of. 


(81) 


For  example,  in  this  notation  the  eddy  viscosity,  A  (Equation  6), 
assumes  the  form 


A  — 

where  « is  the  “dimensionless  eddy  coefficient.” 

The  Mz-tide  is  computed  with  the  reference  values 


Sr  =  30°,  Hr  =  7259.84  m. 


so  that 


y,  =  3/4,  ttr  =  0.021  919  2  sec~\ 
ftir  =  5,  Arr  =  186.309  sec, 

1 80oA/r/' *■  =  1.5°,  J,  =  360°/1.5c  =  240, 


(82) 


(83) 


(84) 


where  JP  is  the  number  of  time  steps  required  to  integrate  through 
one  tidal  period  with  At  =  At,.  Due  to  the  grading  of  the  grid  sys¬ 
tem  (Equation  49),  one  has  almost  everywhere  (n  *  1,  2,  3,  and 
166,  167,  168) 

1  ,  1  (85) 


Similarly,  due  to  the  cutoff  depth  data  (Equations  50,  51),  one  has 
the  limits 


0.0065  <h  <  1. 


(86) 
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Several  other  reference  values  within  the  stability  limits  have  been 
explored.  Particularly  extensive  computations  were  carried  out  with 
Atr  =  248.412  sec  (so  that  Jp  =  180),  but  the  above  reference 
values  appeared  to  yield  the  best  results. 

With  the  simplifications  and  notations  above,  the  coefficients 
(Equations  62-65)  of  the  homogeneous  DOTEs  (Equations  60 
and  61)  become: 


A1  =  A2  =  B1  =  B2  =  0, 

A3  =  rBs  =  PGHT&t/L, 

A*  =  B*  =  2ktr$(0, 

A*  =  A6  =  B6  =  B*  =  htrf/wr, 
A7  =  A8  =  BT  =  r2A5, 

A9  =  A10  =  B8  =  B°  =  B10  =  0, 
Cx=  Taz/L,  C2  =  C3  =  M/L. 


>  (87) 


The  reduced  DOTEs  (with  constant  coefficients)  yield  the  Fourier- 
type  eigensolutions 


and 


UL.n  =  t/odV  [yi(2m  -  2m)AA  +  ya(2n -  1 ) A0], 
VL.n  =  Vod'e*  [yi(2 m  -  *)**  +  y*(2 n  -  1)A0], 


>  (88) 


~  CodV  [yi(2m  —  m)  aa.  4-  yt(2n  —  1)A0].  J 


with  an  arbitrary  wave  vector  (yi,  y2)  and  some  nonzero  amplitude 
vectors  (Uo,  Vo,  Co),  provided  the  eigenvalue  d  satisfies  the  cubic 
characteristic  equation 


(An  dA  io ) 

0 

(1  -n  +  sd)A$i 


(A 25  ~  dA  10 ) 
(l-*  +  *rf)i4» 


Ait 

An 

(1  ~d) 


=  0, 


(89) 


% 
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where  (after  some  algebra) 


and 


with 


A 10  —  1  +  *A*  =  1  + 

Au  —  An  —  A  io  ~  2A4s3, 

AisAai  +  A12A23  — 

_ „  ^  F2  sin*yiA0  +  sin3y2A0  ^ 

0^  ST  ^  rw  7  "i  ^ 


2*1 


\ 


( 


(90) 


(91) 


The  cubic  characteristic  Equation  89  yields  die  three  eigenvalues 

(do,  d  =  di,  da) 

do  —  1  •—  (92) 

and 

dAio  =  A 10 ~  HtT&fAoi  *  2hs  [h®{PAio  —  h^A\i)]U2, 

(93) 

where 

Am=«$  +  ftZT.  (94) 

The  DOTEs  will  be  stable,  provided 

jdfcj  ^  1  for  it  =  0, 1,  2.  (95) 

Under  die  strict  inequality  of  Equation  95,  the  three  eigenvalues 
do,  di,  and  da  define  three  decaying  eigenwaves  represented  by 
Equations  88.  Since  do  is  real,  the  corresponding  eigenwave  is  a 
standing  wave  with  no  phase  shift  if  do  >  0.  The  other  two  eigen¬ 
values  di  and  da  define  a  pair  of  eigenwaves  progressing  in  opposite 
directions  with  the  same  decay  and  dispersion  rates,  provided 


k 


t 
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PA10^h<S^A20u  (96) 

This  condition  holds  true  for  all  0  ^  s  1  (Equation  91 ) ;  i.e.,  for 
all  wave  vectors  (yi,  ys)  if  and  only  if 

0Au>^A<3A2,„.  (97) 

If  this  condition  fails,  then  there  exist  some  short  waves  with  large 
wave  numbers,  yi  and  yi,  which  become  standing  waves  of  different 
decay  rates.  However,  all  sufficiently  long  waves  remain  dispersively 
progressing  and  decaying  at  the  same  rates. 

It  seems  physically  plausible  to  treat  all  long  and  short  waves 
equally  and,  hence,  to  impose  the  conditions  of  Equations  95  and 
97  on  the  free  parameters, «,  r,  *,  and  tt.  Using  Equations  93,  94, 
95,  and  97,  one  finds 

|di|2=  \d2\2=  1  -4A«m2  [tlf- (1  -*)/3t]/Ai0  (98) 


with 


\dx\*  =  |*|*  =  do  fonr  -  1.  (99) 

For  *  <  1,  the  stability  condition  requires 

«*5*(1-T),8r;  (100) 

i.e.,  a  minimum  of  eddy  viscosity  is  necessary  for  stability.  How¬ 
ever,  for*  =  1,  no  minimum  eddy  viscosity  is  required,  which  ex¬ 
plains  the  choice  made  here  (Equation  66)  and  by  Zahel  (1970, 
1973, 1975)  and  Estes  (1975, 1977). 

For  the  chosen  value  *  =  1,  the  stability  condition  is  satisfied  for 
all,  s,  and  all  eigenvalues  do,  di,  and  da,  when 

AhtJto^Aus,  (101) 


i.e.,  when 
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r  =— —  ^  V4  (2  —  (102) 

**fr 

provided  (Equation  97  in  explicit  form)  also 

+  1  -  *)  +  «*?.  ( 103 ) 

The  obviously  increased  stability  limits  imposed  by  both  conditions 
explain  the  choice  of  «  =  1  made  here  for  the  present  tide  model 
(Equation  66)  in  deviation  from  the  value  *  =  0  used  by  Zahel 
and  Estes. 

With  p  =  0.90  (Equation  17)  and  the  possible  values  of  h  =  1 
and,  simultaneously,  y=to  =  4  (Equations  85,  86),  one  finds  from 
Equations  102  and  103  for  *  =  1  and  r  =  1  the  allowable  range 
for  the  dimensionless  eddy  coefficient,  <: 


0<«<0.3(A/  =  A*r).  (104) 

The  same  range  holds  also  for  the  southern  three  colatitude  lines 
n  =  166, 167, 168,  which  violate  the  condition  of  Equation  49,  but 
for  which  the  relative  depth,  h,  falls  sufficiently  below  unity.  The 
upper  limit  on  <  could  be  raised  somewhat  by  considering  the 
simultaneous  values  of  £  to,  and  h  on  each  colatitude  n  =  4,  5, . . . , 
168  separately.  In  order  to  obtain  the  best  possible  tidal  field,  ex¬ 
tensive  trial-and-error  computations  led  to  the  choice 

At  =  Atr  =  186.309  sec  and «  =  0.075.  (105) 

Tim  final  choice  completes  the  detailed  parameter  specifications  of 
the  M2  tide  model. 

At  this  point,  it  may  be  noticed  that  the  stability  requirement  for 
die  DOTEs  restricts  the  possible  amount  of  eddy  dissipation.  As 
is  physically  plausible,  the  finite-differencing  parameters,  *  and 
x,  the  mesh  size,  L,  and  depth,  H,  the  time  step,  At,  and  the  dimen¬ 
sionless  eddy  coefficient,  «,  are  intimately  related  to  each  other. 
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Trial- and-error  computations  are  needed  to  select  those  parameters 
for  best  results.  It  is  particularly  important  to  observe  that  (espe¬ 
cially  for  l )  the  rates  of  decay  (Equations  92,  98,  and  99)  of 
all  eigenwaves  depend  directly  on  the  product  h<.  Acordingly,  for 
fixed  e,  waves  in  deep  (h=>=l )  ocean  basins  decay  faster  than  in 
shallow  (h«  1 )  regions  if  bottom  friction  is  negligible. 

It  is  obviously  this  physically  realistic  phenomenon  that  led  to 
the  introduction  of  the  novel  depth-dependent  eddy  viscosity,  A, 
defined  by  Equations  6  or  82.  For  a  depth-independent  eddy  vis¬ 
cosity,  one  has  he  =  constant,  in  which  case  waves  would  decay  at 
the  same  rate  in  deep  or  shallow  (see  next  section)  oceans,  even 
though  no  bottom  friction  is  present.  Following  Zahel  ( 1970,  1973, 
1975)  and  Estes  (1975,  1977),  the  present  tide  model  also  used  at 
first  a  constant  eddy  viscosity  with  rather  disappointing  results 
caused  by  the  strongly  varying  bathymetry. 

It  is  interesting  to  note  that  for  the  limiting  case  of  Equation  97  ; 
i.e.,  for 


/Mxo-AJMV,  (106) 


Equation  93  assumes  the  simple  form 


dAn  =  Aoi  -  2/St J*  *  2f/8 rj(  1  -  T1)1'*.  ( 107) 

Hence,  for  the  north-east  waves  under  45°  with  wave  numbers 
(Equation  91), 


yi  =  ya  —  y,  s3  —  sin*  yAQ, 


(108) 


and  (*  =  *  =  1 ) 


«v  +  ft* 


(109) 


Hence,  in  this  case,  the  eigenvalues  dx  and  ds  lie  on  the  circle 
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\d~dc\=dr,dc  +  dr=l  (110) 

as  illustrated  in  Figure  2. 

Discussion  erf  the  Tide  Model 

Discrete  Venu*  Continuous  Ocean-Tide  Equation* 

In  the  following  discussion  it  may  be  contended  that  discrete  ocean 
tide  equations  (DOTEs,  Equations  60  and  61 )  reflect  die  physical 
reality  of  ocean  tidal  currents  more  perfectly  and  perceptibly  than 
the  corresponding  continuous  equations  (COTEs,  Equations  29 
and  30).  Although  the  latter  follow  from  the  former  by  a  “formal” 
limit  process,  it  is  neither  technically  feasible  nor  theoretically  de¬ 
sirable  to  seek  convergence  of  the  discrete  solution  to  the  continu- 


Flgnre  2.  Illustration  of  eigenvalues  d*  dlt  and  da  =  <7„  in  cstle 
I d-d'\*4r. 


ous  integral.  In  fact,  in  fluid-flow  problems  of  global  dimensions,  it 
is  not  possible  to  approximate  the  conditions  of  the  continuous  case 
to  any  reasonable  degree.  Even  with  future  computer  technology, 
it  will  not  be  meaningful  to  refine  significantly,  for  instance,  the 
1°  by  1°  grid  system  defined  in  the  section,  “The  l°xi°  Graded 
Grid  System  and  Bathymetry,”  which,  with  its  100-km  mesh  size, 
is  far  from  being  anywhere  near  a  continuous  description.  Any  at¬ 
tempt  to  refine  the  grid  system  would  have  to  be  matched  by  an 
improved  bathymetry  which  requires  worldwide  in  sLu  measure¬ 
ments. 

From  the  theoretical  point  of  view,  it  is  equally  superfluous  to 
seek  an  (-approximation  of  the  continuous  situation,  which  in  fact 
is  only  vaguely  defined  (section,  “The  Navier-Stokcs  Equations  of 
Averaged  Turbulent  Flow.”)  In  laminar  viscous  flow  theory,  such 
notions  as  particle  (point)  velocity  and  pressure,  as  well  as  the 
Navier-Stokes  equations,  are  all  derived  from  physically  sound  dis¬ 
crete  (finite-difference)  definitions  by  a  formal  limit  procedure;  i.e., 
by  simply  assuming  the  existence  of  the  limit  values  (see,  e.g., 
Schlichting,  1968;  Whitaker,  1968).  While  this  assumption  is  well 
justified  in  most  laminar  flow  problems,  it  is  well  known  today  (see, 
e.g.,  Ladyzhenskaya,  1969)  that,  in  general,  even  laminar  motions 
must  be  sought  in  the  class  of  generalized  (distribution)  functions. 
Hence,  velocities,  pressure,  and  their  derivatives  do  not  exist  in  the 
ordinary  sense  (pointwise);  only  their  “functionals”  (effects  such 
as  mass  fluxes,  forces,  momenta)  are  physically  defined. 

The  ambiguity  of  continuous  flow  models  becomes  much  more 
apparent  in  the  critical  laminar  regime.  Experiments  (e.g.,  Busse 
and  Whitehead,  1971)  and  theory  (e.g.,  Schwiderski,  1972)  have 
clearly  established  that,  when  given  characteristic  flow  parameters 
(dimensions,  velocities,  etc.)  exceed  certain  critical  values,  the 
corresponding,  uniquely  existent  laminar  motions  become  unstable 
and  bifurcate  into  laminar  flows  in  infinitely  many  different  shapes. 
The  rUffltffcy*1  laminar  boundary  and  initial  conditions  are  no  longer 
sufficient  to  specify  a  unique  motion.  The  situation  seems  to  be 
governed  by  hysteresis  and  pure  chance  rather  than  by  rigorous 
physical  selection  principles. 
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Critical  laminar  motions  become  still  more  microscopically  un- 
definable  when  die  corresponding  characteristic  flow  parameters 
(as  die  global  dimensions  of  ocean  currents)  exceed  further  super¬ 
critical  points  and  the  motions  go  turbulent  The  statistical  ap¬ 
proach  underlying  the  “time-averaging”  process  to  derive  die  so- 
called  Navier-Stokes  equations  of  mean  turbulent  flow  (Equations 
1  and  2)  is  entirely  formal  and  vague  (see,  e.g.,  Schlichting,  1968). 
For  example,  what  velocity  (particle,  point  etc.)  is  averaged  over 
what  time  interval?  In  this  respect  the  present  periodic  tidal  mo¬ 
tions  are  clearly  die  most  illuminating  of  the  problems  at  hand.  If 
one  averages  (as  usually  meaningful)  over  a  “sufficiently”  long 
time  span  (say,  longer  than  the  tidal  periods),  then  the  averaged 
velocity  should  approach  zero,  which  is  obviously  not  of  interest. 

Evidendy,  turbulent  particle  velocities  manifest  themselves  sta¬ 
tistically  through  their  integrated  (macroscopic)  physical  effects, 
such  as  mass  fluxes.  Hence,  the  proper  mathematical  representation 
of  turbulent  motions  should  be  sought  in  the  class  of  generalized 
functions.  Since  the  product  of  generalized  functions  has  no  mathe¬ 
matical  meaning  (see,  e.g.,  Shilov,  1968),  it  appears  understand¬ 
able  that  there  is  no  way  to  define  the  Reynolds  stress  tensor  of 
turbulent  motion  (section,  “Reynolds  Stresses  and  Eddy  Dissipa¬ 
tion”)  by  a  meaningful  ordinary  or  generalized  function,  because 
it  contains  quadratic  products  of  the  so-called  fluctuating  velocity 
residuals  (see,  e.g.,  Schlichting,  1968).  However,  its  energy-dissi¬ 
pating  (stresslike)  effect  is  physically  quite  apparent  and  must  be 
modeled  in  some  macroscopic  sense. 

To  avoid  all  conceptual  difficulties  of  microscopic  turbulent  mo¬ 
tions,  it  seems  natural  to  fall  bade  to  the  discrete  (macroscopic) 
description  of  laminar  flows  that  leads  formally  to  the  Navier- 
Stokes  equations.  In  fact,  by  proper  “generalized”  interpretation, 
most  motions  used  in  the  laminar  regime  retain  their  physical  mean¬ 
ing  in  the  discrete  turbulent  domain.  For  example,  the  “mean” 
x-vdocity,  u,  of  a  “flow  parcel”  contained  in  a  rectangular  test 
(grid)  cell  of  mesh  lengths  *x,  Ay,  and  Az  (Figure  3a)  at  some 
time,  t,  of  a  time  interval.  At,  is  defined  as  the  mass  flux  (AM*) 
crossing,  say  tt «  central  (x  =  x»)  surface  element  (Ay,  Az)  of  the 
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Ax/2  xf  Ax/2 


Figure  3a.  Illustration  of  average  velocity:  4WJ  =  mass  flux  in  x-dircetion, 
u  =  mean  velocity  in  x-direction,  and  *  =  fluid  density. 

cell  during  die  time  span  At  divided  by  die  fluid  density,  p,  the  area, 
Ay  •  Az,  and  the  time.  At,  so  that 

ii  =  AAfVpAyAxA/.  (Ill) 

This  (generalized)  definition  of  the  mean-flow  velocity,  u,  is 
uniformly  valid  for  laminar  and  turbulent  motions,  since  AM*  is  a 
physically  realistic  (measurable)  quantiy  for  any  p  >  0,  Ax  >  0, 
Ay  >  0,  Az  >  0,  and  At  >  0.  If  u  exists  in  the  limit  Ax-*0,  Ay-»0, 
Az-»0,  and  M-+0,  then  u  becomes  the  ordinary  particle  (or  point) 
velocity.  This  assumption  is  well  justified  in  most  (see  above) 
laminar  motions,  where  neighboring  particles  flow  on  smooth 
neighboring  pathlines.  In  contrast  to  laminar  flow,  the  limit  assump- 
ion  is  certainly  not  justified  in  turbulent  flow,  where  the  fluid  parti¬ 
cles  move  in  an  indeterministically  turbulent,  “fluctuating,”  or 
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“eddying”  way.  In  the  sense  of  generalized  functions,  u  could  be 
defined  in  the  limit  by  prescribing  its  mass  flux  (AM1)  for  every 
nonvanishing  Ax,  Ay,  Az,  and  At. 

Using  similarly  generalized  definitions  of  all  mean  velocity  com¬ 
ponents  and  pressure,  and  retaining  their  finite  differences  in  place 
of  (generally  non-existent)  ordinary  derivatives,  one  arrives  as 
usual  at  the  discrete  Navier-Stokes  equations  (analog  of  Equations 
I  and  2)  without  requiring  the  existence  of  any  limit  values  (see, 
e.g.,  Schlichting,  1968;  Whitaker,  1968).  They  express  simply  the 
physical  laws  of  conservation  of  momentum  and  mass  for  every  test 
(mesh)  cell.  This  is  particularly  tangible  for  the  DOTE  (Equation 
61),  which  conserves  mass  by  balancing  the  excess  of  mass  flux 
into  a  mesh  cell  with  a  corresponding  increment  in  tidal  height. 

In  the  discrete  form,  the  Navier-Stokes  equations  look  formally 
the  same  for  laminar  and  turbulent  motions.  Yet,  when  in  the  lam¬ 
inar  case  velocity  and  pressure  exist  pointwise  in  the  limit,  one  has 
a  unique  microscopic  (continuous)  description  of  the  flow  subject 
only  to  specified  boundary  and  initial  conditions.  Hence,  the  qual¬ 
ity  of  a  solution  of  the  discrete  model  can  be  measured  against  the 
so-called  exact  integral  of  the  continuous  equations.  In  sharp  con¬ 
trast  to  the  laminar  case,  no  unique  microscopic  description  of  tur¬ 
bulent  flow  exists  against  which  the  quality  of  the  discrete  model 
could  be  measured.  Therefore,  die  numerical  analyst  must  seek  an 
optimum  size  of  the  mesh  cell  (Ax,  Ay,  Az)  and  the  time  span  (At) 
in  order  to  model  the  undetermined  microscopic  turbulent  motion 
so  that  their  macroscopic  effects  match  the  expected  or  observed 
features. 

The  strong  dependence  of  a  discrete  turbulent-flow  model  on 
the  size  of  the  mesh  cell  (Ax,  Ay,  Az)  and  die  time  span  (At)  can 
be  assessed,  for  instance,  from  die  definition  of  the  average  ve¬ 
locity,  u,  by  Equation  111.  With  increasing  mesh  area  (Ay  •  Az), 
more  and  more  as  well  as  larger  and  larger  fluctuating  or  eddying 
motions  are  filtered  out  and  remain  unaccounted  for  in  die  average 
value  of  u.  Hence,  the  maximum  mesh  lengths  Ax,  Ay,  and  Az  must 
be  sufficiently  smaller  than  the  smallest  wave  length  one  wishes  to 
resolve.  On  the  other  side,  if  the  area  (Ay  •  Az)  in  Equation  1 1 1  is 
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chosen  smaller  and  smaller,  then  u  becomes  more  and  more  unde¬ 
termined  (fluctuating). 

Similar  arguments  determine  an  optimum  time  step,  At.  In  the 
present  discrete  tide  model,  the  cell  size  was  reasonably  limited  by 
the  available  bathymetric  tables.  The  time  step  At  (Equation  10S) 
was  determined  by  trial-and-error  computations,  so  that  60  time 
points  represent  one-quarter  of  the  Ms-tide  period. 

Another  significant  distinction  between  the  discrete  Navier- 
Stokes  equations  of  laminar  and  turbulent  flow  becomes  apparent 
in  the  average  stress  tensor.  As  is  well  known,  the  turbulent  fluctu¬ 
ations  neglected  in  the  mean  velocity  and  pressure  manifest  them¬ 
selves  as  stresslike  (energy-dissipating)  forces  that  affect  the  mean 
motion.  Unfortunately,  no  exact  and  unique  constitutive  equation 
is  known  today  that  relates  those  turbulent  Reynolds  stresses  to  the 
mean  rate  of  strain  (deformation  of  the  flow  parcel)  determined 
by  the  average  velocity  (see,  e.g.,  Schlichting,  1968;  Whitaker, 
1968).  In  the  Boussinesq  (1877)  substitution  used  in  the  present 
tide  model,  the  mean-turbulent-stress  tensor  is  directly  related  to 
the  average  rate  of  strain  (analog  of  Equations  S).  Hence,  the 
macroscopic  stress  effects  of  the  turbulent  fluctuations  on  the  mean 
velocity  are  assumed  to  follow  a  similar  simple  law  as  the  viscous 
laminar  stresses  in  Newtonian  fluids.  Only  the  coefficient  of  vis¬ 
cosity  is  replaced  by  the  so-called  eddy  viscosity,  which  remains  to 
be  modeled  to  account  for  the  otherwise  neglected  eddying  motions 
in  some  best  sense. 

In  the  absence  of  better  approximations,  it  seems  idle  to  argue 
about  the  physical  justification  of  the  Boussinesq  substitution;  the 
fact  remains  that  it  represents  the  simplest  possible  constitutive 
equation,  including  zero  used  by  some  researchers.  Moreover,  it 
provides  for  considerable  flexibility  to  model  the  microscopically 
undetermined  but  macroscopically  apparent  eddy  dissipation  by 
choosing  suitable  velocity-dependent  or  velocity-independent  eddy 
viscosities  either  uniformly  or  separately  for  all  three  stress  direc¬ 
tions.  Evidently,  if  the  velocity  field  were  known,  a  priori,  then  one 
could  determine  exact  eddy  viscosities,  a  posteriori,  in  many  ways. 
In  this  connection,  it  is  of  interest  to  know  that  the  mean  flow  is 
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quite  insensitive  to  fairly  large  variations  (say  25% )  of  the  eddy 
viscosity.  This  observation  by  Munk  and  Palm6n  ( 1951 )  was  con¬ 
firmed  for  oceanic  tidal  motions  by  the  author’s  extensive  com¬ 
puter  experiments.  It  is  probably  related  to  the  well-known  fact 
that  potential  motions  satisfy  the  complete  Navier-Stokes  equations 
of  laminar  flow  with  any  constant  viscosity.  Above  all,  as  with  any 
other  physical  law,  the  Boussinesq  substitution  has  successfully 
passed  its  crucial  test  in  many  practical  applications  in  hydrody¬ 
namics,  oceanography,  and  meteorology.  The  present  ocean-tide 
model  is  no  exception  (see  Schwiderski,  1979b,  1979c). 

In  order  to  illustrate  die  Boussinesq  substitution  in  the  discrete 
case,  (me  may  consider,  for  example,  the  average  normal  stress  r“ 
produced  by  the  filtered  out  (see  die  remarks  to  Equation  111) 
fluctuating  motions  on  the  surface  (Ay,  Az)  at  x  =  Xi  shown  in 
Figure  3b.  Following  Boussinesq,  (me  has  (see  the  analogous 
Equation  5c) 


FigMr*  3b.  Illustration  of  mean  normal  stress:  u0,  ua  =  average  x- velocities 
at  x  =  x+  xt,  and  rf  =  average  normal  stress. 
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where  u<>  and  ua  are  the  corresponding  mean  velocities  at  x  =  xo  = 
Xi  —  Ax  and  x  =  X2  =  Xi  +  Ax.  Hence,  the  turbulent  stress  *** 
grows  linearly  with  the  rate  of  change  of  average  velocity.  Analo¬ 
gous  to  the  corresponding  laminar  stress,  this  linear  law  appears 
physically  acceptable,  since  the  expected  mean  tidal  velocities  are 
small.  One  concludes  from  Equation  112a  that  a  large  change  of 
mean  flow  produces  a  large  turbulent  stress,  which  plausibly  must 
be  due  to  strong  fluctuating  motions. 

If  one  assumes  a  constant  eddy  viscosity.  A,  in  Equation  112a, 
then  the  analogy  between  turbulent  and  laminar  stress  becomes 
complete.  However,  as  was  pointed  out  above,  the  strength  of  the 
fluctuating  motions  under  consideration  depends  on  the  size  of  the 
surface  area  (Ay  •  az),  which  justifies  the  assumption 

A  — aAyAj  (112b) 


in  Equation  112a,  where  a  may  be  held  constant  or  used  for 
further  modeling.  Similar  arguments  can  be  developed  for  all  six 
turbulent  stress  components  (see  Equations  5)  with  eddy  viscosi¬ 
ties  equivalent  to  Equation  1 1 2b. 

The  novel  eddy- viscosity  law  expressed  by  Equation  112b  is 
obviously  equivalent  to  the  eddy  viscosity  introduced  in  the  present 
tide  model  by  Equation  6a.  It  also  explains  Equation  4a.  which 
specifies  the  bottom-friction  coefficient,  B,  of  the  discrete  tide 
model.  The  need  for  a  mesh-area-dependent  eddy  viscosity  became 
apparent  when  initial  tidal  computations  with  a  constant  value 
failed  to  yield  realistic  results.  It  must  be  emphasized  that  the  mi¬ 
croscopically  indeterministic  nature  of  turbulent  motions  is  not 
completely  removed  in  the  discrete  flow  model.  Its  specific  macro¬ 
scopic  effects  on  the  mean  flow  are  apparent  in  the  required  opti¬ 
mum  choice  of  the  grid  system,  time  step,  differencing  parameters, 
and  most  of  all,  the  eddy-viscosity  law  and  scaling  coefficient.  No 
information  on  the  fine  structure  of  the  turbulence  can  be  expected 
from  such  a  model. 
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Preliminary  M^-Tide  Results 

The  fully  specified  tidal  equations  (DOTEs  60,  61)  have  been 
applied  to  construct  the  principal  semidiurnal  lunar  ocean  tide 
(M2).  The  results  (Schwiderski,  1976)  represented  a  considerable 
improvement  over  earlier  models.  Indeed,  they  demonstrated  the 
feasibility  to  achieve  the  extraordinary  accuracy  of  10  cm  required 
in  many  applications.  For  example.  Goad  and  Douglas  (1977a,  b), 
who  used  the  preliminary  model  to  determine  the  Ms-tide  effects  on 
the  moon’s  orbital  parameters,  found  satisfactory  agreement  with 
other  theories  and  observations. 

Due  to  the  absence  of  an  exact  continuous  ocean-tide  model 
(preceding  section),  the  degree  of  reality  achieved  by  any  approxi¬ 
mate  discrete  model  must  be  measured  against  empirically  known 
features  of  ocean  tides.  Fortunately,  for  the  M2  ocean  tide  a  large 
number  (e.g.,  British.  Admiralty  Tide  Tables,  1977)  of  tidal  ob¬ 
servations  around  the  world  oceans  are  available  for  comparison. 
Table  2  gives  a  statistical  evaluation  of  the  computed  preliminary 
M2-tide  in  comparison  to  observations  made  at  20  island  tide-gauge 
stations  in  each  of  the  Atlantic,  Indian,  and  Pacific  oceans.  The 
agreements  found  are  encouraging. 

Table  2 

Comparison  ol  preliminary  M,  ocean  tide  model  with  empirical  data. 


Amplitude  Phase 

Number  of  Differences  Differences 


Ocean 

Islands 

Mean 

RMS 

Mean 

RMS 

Atlantic 

20 

—0.6  cm 

10.9  cm 

8.8° 

14.4° 

Indian 

20 

+1.9  cm 

15.1  cm 

5.5° 

19.7° 

Pacific 

20 

+7.4  cm 

12.9  cm 

2.8° 

17.4° 

While  the  preliminary  Ms  ocean  tide  model  displayed  good  agree¬ 
ment  with  most  tidal  observations  over  large  ocean  areas,  signifi¬ 
cant  shortcomings  still  persisted,  especially  over  narrow  ocean 
ridges  such  as  the  Aleutian,  Hawaiian,  Marianas,  and  Caribbean 
ridges.  By  comparing  tidal  observations  on  both  sides  of  such  ridges 
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experimental  tidalists  (see,  e.g.,  Harris,  1904;  Bogdanov,  1961; 
Defant,  1961;  Luther  and  Wunsch,  1975)  very  early  discovered 
drastic  distortion  and  retardation  effects  of  those  bottom  barriers  on 
ocean  tides  (see  the  tables  in  Part  H  of  this  paper,  Schwiderski, 
1979b). 

In  order  to  model  local  tide  distortions,  the  described  tide  pro¬ 
gram  will  be  modified  in  Part  II  by  using  a  new  “hydrodynamically 
defined  ocean  bathymetry”  that  reflects  the  barrier  effects  of  narrow 
ridges  and  other  large  irregularities  of  the  ocean  basin  (Schwider¬ 
ski,  1978a).  Moreover,  a  novel  “hydrodynamical  interpolation 
technique”  will  be  introduced  to  incorporate  directly  over  2,000 
empirical  tide  data  into  the  construction  of  the  tidal  charts.  The 
high  numerical  quality  of  the  final  M2  ocean  tide  model  will  be 
displayed  by  maplike  tables. 

A  resume  of  this  paper  was  presented  at  the  International  Sym¬ 
posium  on  Interaction  of  Marine  Geodesy  and  Ocean  Dynamics, 
held  from  October  10  to  13,  1978,  at  the  Rosenstiel  School  of  Ma¬ 
rine  and  Atmospheric  Science ,  Miami.  Florida. 
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